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We discuss the role coarse-grained models play in the investigation of the structure and ther- 
modynamics of bilayer membranes, and we place them in the context of alternative approaches. 
Because they reduce the degrees of freedom and employ simple and soft effective potentials, 
coarse-grained models can provide rather direct insight into collective phenomena in membranes 
on large time and length scales. We present a summary of recent progress in this rapidly evolving 
field, and pay special attention to model development and computational techniques. Applica- 
tions of coarse-grained models to changes of the membrane topology are illustrated with studies 
of membrane fusion utilizing simulations and self-consistent field theory. 
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1. Introduction 

The incredible complexity of biological systems combined with their immediate importance 
makes them the most recent subject for the application of the coarse-grained models of soft 
condensed matter physics [ E 12 EI • While developed earlier for the elucidation of the statics 
and dynamics of melts and solutions of very long and uniform polymers [ El ISII^ 1711511^ ITUIITT] . 
coarse-grained models seem particularly well suited to the study of biological constituents [ 
1121 113L I14L I15L 1161 117j . such as the heteropolymers DNA and RNA, as well as relatively short- 
chain lipids which comprise all biological membranes. 

Systems of polymers and of lipids share many common features, and exhibit universal collective 
phenomena, those which involve many molecules [ I13j . Examples of such phenomena include 
thermodynamic phase transitions, e.g., the main chain transition in lipid bilayers from a fluid, 
liquid crystalline to a gel phase [ 1181 ITU] , and the lateral phase separation which appears to be 
implicated in "raft" formation in the plasma membrane [ 1201 12H I22| I23j , as well as thermally 
activated processes such as vesicle fusion and fission [ US UHl HH UHl UHl 00] j important 
in endocytosis and exocytosis, and electroporation [ |^ EHl ES EEj used in the micro- 
encapsulation of drugs and drug-delivery systems [ EHl EZl- 

Since collective phenomena involve many molecules and entail large length and time scales - 
10-1000 nm and /xs-ms, respectively - details of the structure and dynamics on short, atomistic 
length scales are often irrelevant, and the behavior is dictated by only a small number of key 
properties, e.g., the amphiphilic nature of the molecule. This imparts a large degree of universal- 
ity to the collective phenomena. These terms are borrowed from the theory of critical phenomena 
[ I38j . However the clear separation in length, time and energy scales assumed by this approach, 
is often missing in membrane systems. Thus the universality of collective phenomena, or the 
ability of coarse-grained models to describe collective phenomena, cannot be taken for granted. 
It is important, therefore, to compare the behavior of different experimental realizations among 
each other and with the results of coarse-grained models. 

In the following we shall highlight some recent developments in this active research area 
in which many new models and computational techniques are being developed. We do not 
attempt to provide a comprehensive overview of this rapidly evolving field, but rather try to 
give an introduction both to the basic concepts involved in creating a coarse-grained model, 
and to the simulation techniques specific to membranes and interfaces. We shall emphasize 
the connection to polymer science whenever appropriate. In particular, we will also discuss 
application of field-theoretic techniques to calculate membrane properties. These techniques 
employ very similar coarse-grained models as do the particle-based simulation schemes, and 
they permit the calculation of free energies, and free energy barriers, which are often difficult 
to obtain in computer simulations. An application of coarse-grained models in the context of 
computer simulations and field-theoretic techniques is illustrated by the study of membrane 
fusion, a choice biased by our own research focus on this area. 

Many important applications are not covered by this manuscript. Most notably we do not 
discuss important progress in the study of collective phenomena exhibited by single molecules, 
as in the folding of a protein [ EH ISl 1^ or the processes that ensue when a protein is inserted 
into a membrane [ I42| 143 1 144| I45L I4(ij . or those exhibited by assemblies of a small number of 
molecules, as in the formation and subsequent function of a channel [ 1471 1481 149j . In our view, 
details of the specific molecular architecture are very important for these processes, and they 
lack the type of universality which undergirds the application of coarse-grained models. 

In the next section we place coarse-grained models in the context of atomistic ones that deal 
with molecular details, and of phenomenological Hamiltonians that do not retain the notion of 
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Figure 1. Illustration of different models used to tackle problems in membrane physics, 
(left) snapshot of an atomistic simulation of a 1-stearoyl- 2-docosahexaenoyl-sn-glycero-3- 
phosphocholine (SDPC, 18:0/22:6 PC) lipid bilayer from Ref. [EO]. (middle) tensionless bilayer 
of a coarse-grained model from Ref. [ |^ (right) snapshot of a randomly triangulated surface 
from Ref. [El]. 



individual molecules. We then discuss briefly a selection of simulation and self-consistent field 
techniques utilized for coarse-grained models of membranes. We illustrate the combination of 
computer simulation and field-theoretic approach with the example of membrane fusion. The 
paper closes with an outlook on further exciting, and open, questions in this area. 

2. Atomistic modeling, coarse-grained models and phenomenological Hamiltonians 

Processes in membranes evolve on vastly different scales of time, length and energy. Con- 
sequently a variety of membrane models and computational techniques have been devised to 
investigate specific questions at these different scales. We divide them roughly into atomistic, 
coarse-grained, and phenomenological models as illustrated in Fig. ^ 

2.1. Atomistic molecular dynamics simulations 

Atomistic models, which describe bilayer membrane properties with chemical accuracy, have 
been successfully utilized to investigate the detailed properties of particular membrane systems 
and lipid-protein complexes. Such models have a longstanding tradition. The first Molecular 
Dynamics simulation of lipid bilayers, carried out in the early 1990's were able to sim- 
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ulate only a small patch of a bilayer, one nanometer in extent, over a very short time interval, 
typically 0.1ns. Since then, much effort has been directed to the improvement of simulation algo- 
rithms, (e.g., multiple time step algorithms [ESI)) and to the implementation of the simulation 
code on parallel computers. Consequently atomistic simulations have advanced significantly. 
Today, there are many complete simulation packages that comprise standard routines for Molec- 
ular Dynamics simulations. Among them are NAMD Molecular Dynamics Software (NAMD) [ 
I56j . the Groningen Machine for Chemical Simulation (GROMACS) [ 1571 158j . Groningen Molec- 
ular Simulation (GROMOS) [EHiEni, MDynaMix [|Mj, Assisted Model Building with Energy 
Refinement (AMBER) [EH, NWChem [ESI El, Integrated Model Program, Applied Chemical 
Theory (IMPACT) [ESI; Biochemical and Organic Simulation System (BOSS)/ Monte Carlo for 
Proteins (MCPRO) [ES], DL_POLY [EH EH], Large-scale Atomic/Molecular Massively Parallel 
Simulator (LAMMPS) [ I69j. and Extensible Simulation Package for Research on Soft Matter 
Systems (ESPResSo) [EH- The first four packages in this list are tailored to simulate lipid 
bilayers and proteins in atomistic details, while the latter codes stem from polymer simulations 
and have also been applied to coarse-grained simulations of soft matter. Most of the program 
packages (NAMD, GROMACS, DL_POLY, LAMMPS and ESPResSo) are freely available. A 
nominal fee is charged for the use of GROMOS, while other packages, (e.g., IMPACT), are 
commercial products. Most of the programs run on parallel computing platforms using the Mes- 
sage Passing Interface (MPI). Differences exist in the way the information is distributed among 
the processors (e.g., spatial domain decomposition vs. particle decomposition), the force fields 
that are implemented, and their ease of use and the availability of tutorials. In addition to 
ensuring that the code is error- free and computational efficient, the developers make eff'orts to 
keep it transparent and extendable. Often the data formats allow exchange of information with 
other software packages, such as high level quantum chemistry packages (QM/MM methods), 
or visualization software (e.g.. Visual Molecular Dynamics (VMD) [|^ or PyMol [I72j). 

In atomistic models the attempt is made to describe faithfully the molecular architecture and 
the interactions between components. The interactions include those specifically describing a 
chemical bond (two particle, bond angle potentials and dihedral, torsional, potentials), as well 
as those between atoms not bonded to one another, such as electrostatic and van der Waals 
interactions. The quality of the interaction potentials is crucial for a successful description. In 
lieu of a first-principles calculation, one would like the interactions to fit the results of quantum 
chemical calculations of small fragments of the molecules. Often parameters are additionally 
fitted to experimental observations of various quantities. Many sets of potentials have been 
devised for lipids in an aqueous environment, but there is still a need to refine the interactions 
and devise more accurate models. The parameters of the interactions can be adjusted to standard 
force fields, e.g., CHARMM [17S|, GROMOS96 [TT, OPLS-AA [ESlIZni, or Encad [^77. Often 
there exists the possibility of including customized, tabulated, potentials. 

Atomistic force fields typically include Coulomb interactions. They arise from ionic groups or 
"partial charges" that account for the ability of atomic species to share charges on a common 
bond. For very small systems, the long-ranged Coulomb interactions are handled by Ewald 
summation. Most modern applications, however, employ Particle-Mesh-Ewald techniques [ I78| 
1791 l8Uj which yield a scaling of the order 0{N In N), with the number of charges, N, or fast 
multipole expansions (e.g., IMPACT [ES]) which scale linearly with A^. ESPResSo [ l7Uj . however, 
can additionally deal with Coulomb interactions via a strictly local, field-theoretic, algorithm 
[ 1^ IS21 EHj , and provides routines for simulating systems that are periodic in one or two 
directions. 

Typically, atomistic models are studied by Molecular Dynamics simulations [ El] . The ad- 
vantage of this simulation scheme consists in its rather realistic description of the microscopic 
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dynamics of the constituents. It thereby permits the investigation of kinetic processes, e.g., 
the transport of small molecules across the bilayer, the lateral self-diffusion of lipid molecules 
in the bilayer, the tumbling motion of the lipid tails, or the dynamics of the hydrogen-bond 
network of water at the hydrophilic/hydrophobic interface. The simulation packages often use 
the time-reversible, symplectic Velocity Verlet or Leapfrog algorithms [ 1851 186j and allow for 
multiple time step integration [ 55] . 

Most of the atomistic simulation packages also include a limited selection of methods to cal- 
culate free energies. The most popular techniques comprise thermodynamic integration and 
umbrella-sampling methods [IHZI- Sometimes replica-exchange Monte Carlo, or parallel temper- 
ing methods, are employed [ ES ES EOl EU • 

In the simplest case, the simulation follows the trajectory of the particles by integrating 
Newton's equations of motion in time. The time step is set by the shortest period in the 
system, e.g., that of the stretching mode of a covalent bond. Since the time scale for these 
bond vibrations is orders of magnitude smaller than the time scale of interest, the bond lengths 
are often constrained, thus eliminating the smallest periods and allowing for a larger time step 
of the integrator. Common algorithms that incorporate these constraints into the Molecular 
Dynamics scheme are SHAKE [ 92 , RATTLE [113, and LINCS [El. Moreover the simulation 
packages also allow the possibility of constraining the position of atoms, or distances between 
constituents, a constraint often applied during the relaxation period. 

Newton's equations of motion lead to a microcanonical trajectory. Unfortunately, most nu- 
merical integrations do not conserve the energy on long time scales. Therefore one often couples 
the system to a thermostat, and uses the temperature, T, as a control variable. Analogously, it 
is often desirable to simulate the system at constant pressure, or more usually, surface tension, 
using the Berendsen thermostat [ f0F or the Anderson-Nose-Hoover algorithm [ I96L 19 7| I98| I99j . 
This allows the volume, or area, to fluctuate. Some programs also permit the geometry of the 
simulation box to change in order to equilibrate stresses in crystalline phases via the Parrinello- 
Rahman technique [ llUUj . These simulation methodologies are very similar to what is utilized 
in coarse-grained simulations. 

In addition to interactions between atoms, external forces have also been included in recent 
simulation schemes. In these "steered Molecular Dynamics simulations" one can mimic, for 
example, the action of an AFM cantilever on a biopolymer. This allows one to "push" the 
system over a free energy barrier. If this is repeated often enough, Jarzynski's equation [ llUll 
ll()2( ll(K-i| 11041 ll()5j can be used to calculate free energies by integrating the corresponding 
Boltzmann factor. 

While atomistic simulations in the early 1990s could simulate only an extremely small patch 
of a pure lipid bilayer membrane on the order of a nanometer, [ 1531 I54j . current Molecular 
Dynamics simulations study membrane patches of a few tens of nanometers over time scales 
of a few tens of nanoseconds. The time scale of observation has increased by two orders of 
magnitude over the earliest attempts. These simulations provide valuable information about 
the molecular organization and dynamics of the lipids and of the water, both in the bilayer and 
at the hydrophobic/hydrophilic interface [ 106} I107| I108| \T09\ . For example, recent atomistic 
simulations revealed the molecular structure of the ripple phase in phosphatidylcholine bilayers 
[ IllOj . However, calculation of these basic properties (e.g., free energy differences between 
different morphologies) still pose challenges to computer simulations because of the large time 
scales which are introduced by the interactions between head groups due to charge pairing or 
water bridging. 

Another important aspect of atomistic simulations is that they are able to correlate the de- 
tailed lipid architecture (e.g., particular head group structure, length, and degree of saturation 
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of the tails) with the physical properties of the bilayer (thickness, orientation of segments, liquid- 
crystalline ordering, and elastic moduli). Such physical properties can be obtained, for example, 
by increasing sufficiently the length and time scale of the simulation, so that one can observe 
the fluctuations of the bilayer membrane. From an analysis of the undulations and peristaltic 
fluctuations of the local thickness, one obtains the bending rigidity and area compressibility of 
the membrane, as well as the pertinent relaxation times [ 11111 rrT2| . For instance, atomistic 
simulations of monoolein, which is an extremely simple lipid with only one mono-unsaturated 
tail of eighteen carbon atoms and a very small headgroup consisting solely of two hydroxyls, 
reveal that fluctuations on the length scale of 20 nm have relaxation times of more than 5 ns. 
This time scale is expected to be longer in the usual double-tailed phospholipids [ I112j . 

The combination of sophisticated simulation codes and powerful parallel computers has per- 
mitted atomistic simulations to investigate some collective phenomena in bilayer membranes. 
The transformation from an inverted cubic phase to an inverted hexagonal phase in the monoolein 
system has been studied [ I113j . The transition, as a function of areal density, between liquid- 
condensed and liquid-expanded phases of a monolayer of DPPC at an air-water interface has 
been studied, as well as the monolayer's eventual rupture [ I114j . Recent applications can even 
cope with the protracted time scales associated with the spontaneous self-assembly of lipids into 
a bilayer [ I115j , or the formation of small vesicles [ I116j . 

Another thrust of applications is the study of membranes that consist of lipid mixtures, or 
membranes that have molecules adsorbed onto them. The atomistic description provides a 
detailed view of the role that small inclusions, like cholesterol, or adsorbents, like sugars or 
polymers, play in modifying the structure of the bilayer. In the following, we give several 
examples. 

The addition of a surfactant [ I117j or cholesterol [ I11S| I119| I12()j to a lipid bilayer tends 
to increase the ordering of the lipid tails, which causes them to lengthen. Because the liquid- 
like interior of the bilayer is highly incompressible, this results in a decrease of the area per 
head group. In the case of high concentrations of cholesterol [ I119j . the decrease of the area 
per head group is accompanied by an increase of the bending modulus and a decrease in the 
lateral self-diffusion coefficient of the lipids. The change of the lipid packing also affects the 
distribution of voids. They become rarer, and those that remain elongate along the bilayer 
normal as the concentration of cholesterol increases [ ll2Uj . Addition of salicylate to a lipid 
bilayer also decreases the area per head group. However in this case, the mechanical properties 
of the bilayer are hardly affected [ I121j . The effect of halothane, an anesthetic, is quite different. 
This small molecule preferentially segregates to the upper part of the lipid acyl chains, increases 
the area per head group, and decreases the lipid chain order [ I122j . 

Large sugar molecules do not penetrate into the hydrophobic interior of the bilayer, but 
do impact the hydrogen bonding at the interface between the head groups and the water. An 
interesting example is provided by trehalose, a disaccharide, which is found in animals capable of 
enduring cold temperatures or dry environments. Experiments indicate that it prevents leakage 
and fusion during drying and freeze-drying, a property which has been exploited for practical 
apphcations [ 123, 124,. Atomistic simulations [USSl lEEl IlSZl lEHl Il29,, 130^ show that the 
area per head group remains unaffected. In addition, the total number of hydrogen bonds of 
the lipid heads is conserved. However, hydrogen bonds with trehalose now replace some of the 
bonds which had been made with water. A single trehalose molecule can interact with multiple 
lipids simultaneously. This result suggests that, at sufficiently high concentrations, disaccharides 
might serve as an effective replacement for water. 

The largest-scale simulations carried out on the atomistic level are able to study lipid-protein, 
or lipid-DNA, interactions [ I131j . and to investigate channels [ I47 ( 1132] through a bilayer lipid 
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membrane. The added complexity brought about by incorporating proteins into the membrane 
poses huge chahenges to both the simulation techniques and computational requirements due 
to the large number of additional interactions which have to be accurately described. Some 
examples of the systems studied are as follows. 

A protein's conformations when it is inserted into the membrane, and the distortion of the lipid 
bilayer in its vicinity, can be studied by atomistic simulation. The protein's interactions with 
the lipids are both strong, compared to the thermal energy scale, ksT, and specific. They are 
difficult to simplify, with the result that the details of the complex architecture on the molecular 
level have to be considered for a quantitative description. Proteins can create interactions within 
the bilayer due to the strain field generated by a mismatch of its hydrophobic region with that of 
the bilayer in which it is embedded. Large molecules like DNA [ ll^-ilj and proteins [l47 |ll,S2|ll,S^-ij 
also give rise to significant interactions outside of the bilayer. 

Often proteins are not isolated in the lipid membrane, but aggregate to structures such as 
pores, or ion channels [ 1461 1471 148j . The detailed structure of these channels has attracted great 
interest in understanding how they function to let some ions pass while stopping others. It has 
been possible to study the permeation of water through an acquaporin pore [ HTl 1132] . These 
simulations reveal the motion of a single water molecule as it passes through the channel. The 
trajectories provide insights into the specificity mechanism by which the channel allows water, 
but not ions, to pass. Recently, the permeability for water and ions of the a-hemolysin/lipid 
bilayer complex has been studied by large-scale computer simulations involving 300 000 atoms [ 
I133j . The application of external electrical fields permitted the ion permeability to be obtained 
as a function of bias voltage, and the selectivity of a-hemolysin to chlorine ions to be elucidated. 

The bilayer structure in almost all of these atomistic simulations has to be pre-assembled 
because the time scale of self-assembly from a homogeneous mixture of lipids and water typically 
far exceeds the simulation time scale. (For an exception, see Ref. [ I115| .) This leaves unanswered 
the question of the thermodynamic stability of the pre-assembled membrane. Even though the 
atomistic potentials are parameterized from the interaction of atoms, the manner in which 
these potentials determine the global stability of the different lipid morphologies is subtle, and 
unknown. Furthermore, the transitions between these lipid morphologies, and the formation of 
out-of-plane structures as occurs in budding, are beyond the scope of atomistic modeling. 

2.2. Coarse-grained models 

2.2.1. Why are coarse-grained models useful? 

While atomistic simulations provide valuable, detailed, information about the local structural 
properties of lipid membranes, they cannot access the time and length scales involved in collective 
membrane phenomena, which are milliseconds and micrometers respectively. One strategy to 
overcome this difficulty is to eliminate some of the detail by coarse-graining the description. 
Coarse-grained models do not attempt to describe the large scale phenomena by starting from 
the smallest atomic length scale, but rather by lumping a small number of atoms into an effective 
particle [ [El [HI HSl CMl Il35, 13^1 IIHZi 1138, iM CHI ESI EHl EH. These particles 
then interact via coarse-grained, simplified, interactions, ones which typically do not include 
electrostatic and torsional potentials, for example. The reduced number of degrees of freedom, 
and the softer interactions on the coarsened scale lead to a significant computational speed-up 
with the consequence that larger systems and longer time scales are accessible. This makes 
possible the study of collective phenomena in membranes, a study not possible via ab-initio 
methods now, or in the foreseeable future. However the loss of chemical detail limits some of 
the predictive power of coarse-grained models. 

The objectives of mesoscopic modeling are twofold: on the one hand, to help to identify those 
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interactions which are necessary to bring about collective phenomena on a mesoscopic scale, 
such as self-assembly, and on the other to elucidate the universal behavior on the mesoscopic 
scale itself. These include the role of thermal fluctuations, or the existence of phase transitions 
between self- assembled morphologies. Mesoscopic models are also an ideal testing ground for 
phenomenological concepts. 

The obvious question which presents itself is how the coarse-graining is to be achieved. What 
are the relevant degrees of freedom and interactions to be retained at the coarse-grained scale in 
order to incorporate the essential physics of the system? This is a fundamental question of any 
model-building procedure which must be addressed when one abandons ab-initio calculations. 

One can respond that, due to the experimentally observed universality of self-assembled struc- 
tures, any coarse-grained model that includes the relevant interactions will capture the qualita- 
tive features. Consequently one should use the simplest possible model in order to take maximum 
advantage of the computational benefits of coarse-graining. This is the strategy of minimal mod- 
els which were the first to study self-assembly. They are still very popular. The question that 
remains is just what are the "relevant" interactions necessary to bring about the collective phe- 
nomena observed in experimental systems. Within the framework of minimal models, one can 
start with a simple model and successively augment it with additional interactions until the 
phenomenon of interest is captured. While this method appears to be rather crude, it does 
provide insight into which interactions, on the length scales of a few atomic units, are necessary 
to bring about collective behavior in membranes. It also contributes to identifying those mech- 
anisms that underly the phenomena and the degree of universality. Alternatively one can try 
to "derive" coarse-grained models from a specific atomistic system, a procedure which is termed 
"systematic coarse-graining". We shall discuss both techniques in turn. 

2.2.2. Minimal models 

The idea of successively eliminating degrees of freedom from a specific mixture of lipid and 
water to "derive" a coarse-grained model is a beautiful and potentially powerful concept. This 
concept of coarse-grained models has a long-standing tradition in polymer physics and dur- 
ing the last three years much progress has been made in the area of biological and synthetic 
membranes. Unfortunately, the coarse-graining procedure is often impractical to implement ex- 
plicitly. Notable exceptions are dilute and semi-dilute polymer solutions for which the concept of 
coarse-graining can be formulated in terms of a consistent theory, one which has been extensively 
exploited [EHHT]. 

The configurations of long, flexible, linear polymers in dilute or semi-dilute solutions are 
characterized by a self-similar, fractal structure. This self-similarity extends from the structure 
of a few monomeric repeat units to the size of the entire molecule, which is comprised of hundreds 
or thousands of monomers. For long chain molecules, there is a clear separation between the 
structure on the length scale of a monomeric unit, which strongly depends on the chemical 
structure and details of the interactions on the atomic scale, and the mesoscopic structure of the 
entire molecule. Clearly the chain dimensions depend on the chemical identity of the monomeric 
units in a very subtle manner, but for the description of the large scale properties a single, 
coarse-grained, parameter, the end-to-end distance, suffices. The background of this statement 
is the observation of de Gennes, in 1972, that the behavior of a long, self-avoiding, walk is 
intimately related to the properties of a critical point in a n-component field theory in the limit 
n ^ [ I145j . This opened the way for the use of the machinery of the Renormalization Group 
for the description of polymer solutions, and placed the heuristic observation of the universality 
of the behavior of long chain molecules within a rigorous theoretical framework. The inverse 
chain length plays the role of the distance to the critical point. The behavior at the critical 
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point is universal, i.e., it does not depend on the microscopic interactions but only on a few, 
relevant, features that characterize a universality class. For polymer solutions the relevant 
properties are the connectivity of the monomeric units along the backbone and the excluded 
volume interaction between monomeric units. By virtue of universality, any model characterized 
by these two properties will capture the behavior of polymer solutions in the limit of long chain 
lengths. The theory has provided detailed insights into the large scale chain conformations in 
dilute and semi-dilute solutions, and has been utilized to describe quantitatively the screening 
of the excluded volume interactions, and the cross-over from dilute to semi-dilute solutions [ 

ElElIZl. 

Biological systems do not exhibit the sort of scale invariance that lies at the heart of the Renor- 
malization Group approach to polymer systems. In particular, there is no parameter, like the 
chain length, that tunes the separation between the microscopic scale of the atomic interactions 
and the mesoscopic structure. Another practical complication is that, in contrast to polymer 
systems in which one considers systems of very few components and with simple interactions be- 
tween them, biological systems are composed of many different, complex, structural units which 
are connected by means of several different interactions. As a consequence, the development of 
coarse-grained models for membranes is more an art than a science. It is often guided by phys- 
ical intuition, computational constraints, and a large degree of trial-and-error. The underlying 
assumption is that, just as in polymer solutions, the qualitative behavior of the membrane de- 
pends only on a few coarse-grained parameters that characterize the relevant interactions of the 
mesoscopic model. This assumption is not justified by a rigorous formal theory. Consequently 
it is a priori unknown what the relevant interactions are that have to be incorporated in order 
for a coarse-grained model to faithfully capture the behavior on mesoscopic length scales. The 
answer to this crucial question depends on the specific system, and on the properties in which 
one is interested. For example, the experimental fact that systems which differ chemically a 
great deal, such as biologically relevant lipids in aqueous solution and amphiphilic water-soluble 
polymers, do exhibit a great number of common morphologies implies that the existence of these 
morphologies can be traced back to a small number of simple interactions. 

A key ingredient is the connectivity of two strongly repelling entities within a single molecule. 
Since these two parts are joined together they cannot separate and form macroscopic domains, 
but organize instead into supermolecular structures so as to minimize the unfavorable contacts 
between their parts. The particular physical mechanisms that cause the repulsion between the 
two entities appear to be less important. 

Another significant experimental observation is the correlation that exists between the gross 
amphiphilic architecture of the components of the system and the system's phase behavior. 
Not only does the size of the molecule set the scale of the self-assembled structures, such as 
the bilayer thickness, but also its "architecture", the relative volumes of the two antagonistic 
molecular parts, can be directly correlated with the various morphologies. This correlation has 
been stressed by Israelachvili [ 11461 1147j . Lipids in which the head and tail groups are of similar 
volumes tend to form bilayers. If the lipid tails volume is enlarged or the headgroup reduced (e.g. 
by replacing a phosphatidylcholine with a phosphatidylethanolamine) then the lipids are said 
to be "cone-shaped" and they tend to form inverted hexagonal phases. In this phase, the lipids 
assemble into tubes with the heads directed inward and the tails outward, and the tubes form 
a periodic hexagonal lattice. This concept also carries over to AB diblock copolymers which 
consist of two blocks composed of Na and Nb, monomeric units. In this case, the fraction, 
/ = Na/{Na + Nb), of one block is employed to parameterize the molecular architecture, and 
it also correlates with the observed phase behavior. From these observations, one can conclude 
that it is crucial to conserve the basic geometry of the molecules during the mapping onto a 



Coarse-grained models for biological and synthetic membranes 



11 



coarse-grained model. 

Notwithstanding these important universal aspects, the details of the molecular architecture, 
interactions and the mechanisms of self-assembly, do vary from system to system. In block 
copolymers, for example, the geometrical conformations of polymers strongly fluctuate and, 
therefore, the average geometrical shape of a diblock copolymer is strongly affected by its en- 
vironment. The balance between the repulsive interaction energy of the two components and 
the conformational entropy that describes the change of available molecular conformations dic- 
tates the self- assembled morphology. In lipid systems, however, the molecules are short and 
rigid. The reduced number of molecular conformations severely limits their ability to alter their 
average geometric shape to adapt to the environment. Thus, the concept of packing rigid, wedge- 
shaped, objects is useful, i.e. the "shape" of the molecules does not depend significantly on the 
environment. A mismatch between the molecular geometry and packing constraints cannot be 
completely accommodated by changes of molecular orientation, so that this mismatch also alters 
the local fluid-like packing. It is this interplay between universal and specific aspects that make 
the development of coarse-grained models in biological systems a challenging one. 

The amphiphilic nature of the molecules and the important molecular geometry are char- 
acteristic of self- assembling systems. These two relevant properties must be captured by a 
coarse-grained model. They differ in detail as to how these properties are incorporated, and 
they have been augmented by additional interactions to provide a more detailed description of 
specific systems. 

One of the simplest self-assembling system is that of oil and water and amphiphile, and many 
simple lattice models were introduced to study it [ ll,S4| I148| I149j . Larson was one of the first 
to ask how some of the simplest specific features of the amphiphile, such as the presence of a 
multi-atom hydrophobic tail and the relative volume of head and tail units, would affect the 
phase structure. While water and oil were represented by a single site on a lattice, amphiphiles 
were modeled as a linear string of nearest or next-nearest sites. Interactions between hydrophilic 
units, water or lipid heads, and hydrophobic units, oil or lipid tails, were described by square- 
well potentials that extended over the nearest and next-nearest neighbors. Like units attracted 
each other while unlike units repelled each other. Monte Carlo simulations of this model yielded 
information about possible phase diagrams of ternary water, oil, and amphiphile solutions [ I134| 
I15UI I151j . This simple lattice model was even able to reproduce the complex gyroid morphology 
that has been observed both in amphiphile solutions and block copolymers. Not surprisingly 
the Larson model shares many features with simple lattice models that have been utilized to 
study the universal characteristics of polymer solutions and melts. In the latter context, simple 
lattice models have been employed to problems ranging from the scaling properties of isolated 
chains in good solvent [ I152j , the equation of state of solutions and mixtures [ 11531 1154j , to the 
ordering of diblock copolymers [ I155| I156j . 

To study further how microscopic details affect macroscopic behavior, one must flesh out these 
schematic models by various structural details. Unfortunately, it is difficult to capture details 
of the geometric shape of the amphiphiles in simple lattice models. The restricted number of 
angles between bonds that connect neighboring amphiphilic units makes very difficult a realistic 
description of the rather stiff tails. Further, in lattice models, the head and tail segments 
invariably occupy identical volumes. Some of these difficulties can be overcome by complex 
lattice models, such as in the bond fiuctuation model of Carmesin and Kremer [ 11381 [TF7| I158j . 
In this model, each segment occupies the eight corners of a unit cell of a simple cubic lattice. 
Monomers along the amphiphile are connected by one of 108 bond vectors that are allowed to 
take lengths, 2, \/5, \/6, 3 or \/lO in units of the lattice spacing. This set of bond vectors is chosen 
such that the excluded volume constraint guarantees that bonds do not cross in the course of 
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local random monomer displacements by one unit in one of the lattice axis. Thus, effects due to 
entanglements are captured. The large number of bond vectors and the extended shape of the 
monomers yields a rather good description of continuum space. For instance, the eighty-seven 
different bond angles permit a rather realistic description of stiffness. Artifacts due to lattice 
discretization are strongly reduced, yet the computational advantages of a lattice model (e.g., 
early rejection of trial moves) are retained [HSni EMI- Mor cover, sophisticated Monte Carlo 
simulation techniques have been implemented for lattice models that allow for a very efficient 
relaxation of the molecular conformations and the calculation of free energies. The model can 
be quantitatively mapped onto the standard Gaussian chain model that is used in self-consistent 
field (SCF) calculations (cf. Sec. 12. This allows for a computationally efficient way to explore 
a wide parameter range as well as to calculate corresponding free energies. Amphiphiles have 
been modeled as flexible chains consisting of a hydrophilic and a hydrophobic block. The solvent 
can be described by a homopolymer chain that consists of hydrophilic segments. Like segments 
attract each other via a square well potential that extends over the nearest fifty-four lattice 
sites, while hydrophilic and hydrophobic segments within this range of interaction repel each 
other. The strength of the interaction between the segments controls the free energy of the 
hydrophilic/hydrophobic interface, where as the relative length of the hydrophilic block, /, tunes 
the spontaneous curvature of a monolayer. The model has been successfully employed to study 
self-assembly in diblock copolymers and their mixtures with homopolymers [ I159j . and pore 
formation of bilayers under tension [ ll6Uj . The bending rigidity of a monolayer and tension of a 
bilayer have been measured via the spectrum of interface fluctuations and bilayer undulations [ 
I159( 116(1] . and the fusion of membranes has also been studied [ I16H 1^ within this framework. 

Although this lattice model includes only the bare essentials necessary to bring about self- 
assembly, it is sufficient to describe its universal properties. A mapping of length scale between 
lattice model and experimental realizations can be established by comparing an experimental 
bilayer thickness in nanometers with the bilayer thickness of the model expressed in lattice con- 
stants. Similarly the model's energy scale can be deduced by comparing the experimental and 
calculated free energy of the hydrophilic/hydrophobic interface. Additional characteristics, such 
as the bending rigidity, or the area compressibility modulus, then can be combined in dimen- 
sionless ratios. A comparison of such dimensionless ratios between liposomes, polymersomes, 
and the bond fluctuation model is presented in Table ^ One observes that these mesoscopic 
characteristics do not strongly differ between membranes formed by long amphiphilic diblock 
copolymers and biological lipids in aqueous solution and that the lattice model is able to repro- 
duce the order of magnitude estimate of these properties. Therefore, this table quantifies the 
universality of amphiphilic systems and justifies the use of highly simplified models [lll-ij. 

An alternative procedure to include molecular detail is to use off-lattice models. Clearly these 
models allow for much flexibility in describing the molecular geometry and they can be studied 
by Molecular Dynamics. A generic off-lattice model has been utilized by Smit and co-workers 
[ I168j to elucidate micelle formation. Water and oil molecules are modeled by Lennard- Jones 
particles while the amphiphile is represented by a collection of particles bonded together via 
harmonic springs. The hydrophobic beads form a linear chain tail, while the hydrophilic head 
beads are all bonded to a single, central bead which, in turn, is attached to the tail. This mimics 
the bulkiness of the lipid head. 

Goetz and Lipowsky [ ll.'Sfij employed a model in which the amphiphiles are comprised of a 
single head bead and four tail segments. Water is modeled by a single bead. The interactions 
between the like beads (head-head, water-water, head-water and tail-tail) are of Lennard-Jones 
type with a cut-off at Vc = 2.5a. The energy, e, and range, a, of the Lennard-Jones potential 
set the scales. The hydrophobic interaction between water and tail or head and tail is a purely 
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Table 1 



Structural and elastic properties of bilayer membranes formed by amphiphilic diblock copoly- 
mers, biological lipids and a coarse-grained lattice model, dc - thickness of membrane hydropho- 
bic core in the tensionless state, / - hydrophilic fraction, Co - monolayer spontaneous curvature, 
AA/Ao - bilayer area expansion (critical value for the experimental systems, and the actual 
strain used in simulations), Ka - bilayer area compressibility modulus, - monolayer bending 
modulus, 7int - hydrophilic/hydrophobic interface tension (oil/water tension of 50pN/nm for 
the experimental systems, and A/B homopolymer tension for the simulations). Data on E07 
polymersomes is taken from [ I162j : and on lipids from (a): [ I163j . (b): [ I164j . (c): [ I165 j. and 
(d): [ I166j (see also http://aqueous.labs.brocku.ca/lipid/). Values of dc, Co and Ka for DOPE 
were obtained by linear extrapolation from the results on D0PE/D0PC(3:1) mixture. Values 
of Kb, 7int) and Co for the simulated model were calculated by us using the method of [ I167j . 
From Ref. [I^. 



repulsive soft potential, 

Vsc{r) = 4e(asc/r)^ (1) 

with (Trep = 1.05cj, that is cut-off at rc = 2.5a. The potentials are truncated and shifted so that 
both the potential and the force are continuous at the cut-off: 

ytrunc(r) = I ^^"^ " ^^'"^ " " " - ''^ (2) 



Kj for r > Tc. 

Beads along the amphiphile are bonded together via harmonic springs 

Hond(r) = A;bond(|r| - af. (3) 

The rather large value /^bondO'^/e = 5000 is chosen to constrain the average bond length to a 
value very close to a. Additionally a bending potential of the form 

Vhend = fcbend(l " COs{9)), (4) 

where 9 denotes the bond angle, is included. By increasing the bending stiffness A^bend ^ 5e, one 
can change the conformations from those typical of a very flexible molecule to those characteristic 
of a rigid one. 

A further step in the coarse-graining procedure is to eliminate the solvent particles while 
preserving their effects implicitly. Since the two-dimensional membrane is embedded in a three- 
dimensional volume filled with solvent, the number of solvent particles increases much faster 
than the number of amphiphiles as one studies ever larger systems sizes. Yet the role of the 
solvent often is only to stabilize the bilayer membrane whose properties are the focus of attention. 
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Therefore, if the solvent could be eliminated, it would result in a huge computational speed- 
up. Typically, one can assume that the amphiphile-solvent mixture is incompressible on a 
coarse scale. Then the system configuration is completely described by the configuration of the 
amphiphiles, and the interaction between solvent and amphiphiles can be integrated out giving 
rise to an effective interaction between the amphiphilic units. 

The resultant implicit solvent models have enjoyed long-standing popularity in polymer physics 
where the behavior of polymers in solvents of different qualities often is described by polymers 
in vacuum with effective interactions between the monomeric units. Attractive interactions cor- 
respond to a bad solvent and result in a collapse of the polymer, while repulsive interactions 
correspond to a good solvent because the isolated polymer adopts a swollen, self- avoiding- walk 
like shape [ElEHl- While there exists a formal one-to-one correspondence between the thermo- 
dynamic properties of an incompressible polymer-solvent mixture and the corresponding com- 
pressible polymer model with effective interactions, these effective interactions might not be well 
represented by density-independent pair potentials. For instance, by replacing the repulsion be- 
tween solvent and polymer by an effective attraction between the polymer segments, one might 
observe a much higher local polymer density than in the original incompressible mixture where 
the maximal value of the local polymer density is limited by the incompressibility constraint. 
The differences between incompressible mixtures and effective compressible systems comprised 
of only amphiphiles are even more pronounced when one regards dynamical properties because 
(i) the density variations in the implicit solvent model results in variations in the local mobil- 
ity of the amphiphilic units that are absent in the original incompressible system and (ii) the 
solvent carries momentum, and the concomitant hydrodynamic flow might promote coopera- 
tive re-arrangements. These considerations illustrate that integrating out the solvent degrees of 
freedom, though conceptually straightforward, does involve some practical subtleties. 

Initial attempts to construct solvent-free membrane models using simple pairwise interactions 
were rather unsuccessful. The model of Drouffe et al. [ llTOj represented amphiphiles by single 
beads interacting via a spherical hard-core repulsion and an orientation-dependent short-ranged 
attraction. They found that increasing the attraction between the lipid tails resulted in the 
formation of membranes. These membranes consisted of a single layer of particles and the mem- 
branes were crystalline (gel), i.e. the lipids laterally condensed onto a triangular lattice. This 
solid phase was characterized by the pronounced reduction of lateral lipid diffusion. When the 
temperature was raised the membrane did not form a fluid membrane, but simply disassembled. 
To overcome this difficulty, a many-body interaction was introduced to mimic the hydrophobic 
effect and to stabilize a fluid membrane. Additionally these interactions limited the number of 
neighbors and thereby suppressed three-dimensional aggregation in favor of sheet-like structures. 
These multi-body, or density-dependent, interactions made the calculation of thermodynamic 
quantities rather subtle (see below). 

Noguchi and Takasu [ I171j modeled the amphiphiles by rigid rods comprised of three inter- 
action centers, a head and two tail beads. These beads interact via a rotationally symmetric 
potential but the multi-body character of the attraction of the hydrophobic tail beads is used 
to stabilize the membrane. Particles repel each other via a soft core potential which defines the 
energy scale, e, and the monomeric length scale, a. The potential is of the form 



and it is truncated and shifted at a cut-off 1.3a. The multi-body potential takes the form 




(5) 




for p < p* — 1 

for p* - 1 < p < /J* 

for p* < p 



(6) 
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with parameters p* = 10 and c = 4.75 for the tail bead nearest to the head and p* = 14 and 
c = 6.75 for the tail bead at the end. The smoothed density, p, quantifies the local number of 
hydrophobic particles in a small sphere around the reference particle at position, r, according to 

{1 for r < 1.6(7, 

expi20(./Li.9)j+i for 1.6a < r < 2.2a, (7) 
for 2.2cj < r. 

where the sum over r' includes all hydrophobic segments on other amphiphiles. At small 
smoothed densities, p < p* — 1, the multi-body potential is linear in the density and, thus, 
represents a pairwise attraction between neighboring hydrophobic beads on different molecules. 
At higher densities, the attractive strength levels off and adopts a constant value independent of 
the local density. This feature avoids the collapse of the hydrophobic tails into extremely dense 
structures and thus prevents crystallization. In contrast to the previous model of Drouffe and 
co-workers [ ll7Uj the membranes in Noguchi's model are bilayer, i.e., they are comprised to two 
layers of amphiphilic molecules. 

Wang and Prenkel [ I172j described another variant of solvent-free models with multi-body 
interactions, where amphiphiles were modeled as flexible chains consisting of three coarse-grained 
beads. A bending potential along the backbone was utilized to tune the molecular flexibility. 
They employed a qualitatively similar density dependence of the multi-body term, but used an 
anisotropic weighting function to construct the smoothed density, p. 

The first solvent-free model that resulted in the formation of a fluid bilayer from particles 
that interact via simple pairwise isotropic interactions was devised by Farago [ I173j . In this 
model, amphiphiles consist of rigid, linear trimers comprising one head-bead and two tail-beads. 
The interactions were tuned by a rather lengthy "trial and error" process to make the attraction 
between molecules sufficiently strong to support the stability of the membrane, but still weak 
enough so that the membrane would not crystallize. They are as follows. Let (1) denote the 
hydrophilic head bead and (2) and (3) the hydrophobic beads along the rigid amphiphile that 
are spaced at a distance, a. Beads of the same type interact via a Lennard-Jones potential 



Vnir) = 4e, 



r 



(8) 



Interactions between the head and the first hydrophobic bead are repulsive 



Vu{r) = Aeu(^^y , (9) 
and the repulsion between the head and the end tail bead is even harsher 

/ \ 18 

^i3(r) = 4ei3 (^^j . (10) 
The attraction between different hydrophobic beads, however, has a broad attractive minimum 



V23{r) = 4e23 



2 

<723 \ Cr23 



(11) 



All potentials are truncated and shifted at a cut-off distance 2.5o"33. The potential parameters 
are detailed in Tabled 

Deserno and co-workers [ I174| I175| argued that an increase of the range of the interaction 
is crucial for stabilizing fluid bilayers. They represented amphiphiles as flexible trimers. All 
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Table 2 



Parameters of the interaction potentials in Equations (|8))- (|11|) of Farago's solvent-free model. 
From Ref. [ 117^] . 



beads repel each other via a Lennard-Jones potential of the type of Eq. Q which is truncated 
and shifted at the minimum, rmin = v^<7m, resulting in a purely repulsive potential. The size 
of the tails defines the length scale, (T33 = (T22 while the heads are smaller, an = 0.95(733 and 
the repulsive interactions between head and tails are non-additive, ai2 = (T13 = 0.95(T33. In 
addition to this purely repulsive interaction, hydrophobic tail beads interact with each other via 
an attraction with tunable range, Wc (see inset of Fig. [^I: 

{-€ for r < Vmin, 

_e cos2 (Ml^^) for r^m <r< r^^^ + Wc, (12) 
O for rmin + Wc <r. 



The beads are linked together via a FENE potential 
1 2 

Vbond = -2^bond'^oln 1 



(13) 



To. 

with /cbond = 30e/cr33 and maximal bond length, ro = 1.50-33. There is no bond angle potential, 
but the flexibility is tuned by applying a harmonic spring between the head and the last tail 
bead 

Vhendins) = ^^bcnd (n3 " "io^sf . (14) 

with /cbcnd = lOe/dlg. 

Figure 121 presents the phase diagram at zero lateral tension as a function of rescaled temper- 
ature and the range of the attractive interaction, Wc- If the range of the attraction, Wc, is small 
compared to the effective hard core diameter, (T33, the membrane assembles into a solid sheet 
upon cooling. Only if the range of the attraction is sufficiently large does one encounter two 
transitions upon cooling. As the system is cooled, the amphiphiles first assemble into a fluid 
membrane. Upon further cooling, the membrane crystallizes. In this solid phase the lateral 
mobility of the amphiphiles is strongly reduced. The temperature interval in which the fluid 
membrane is stable increases with the range of the attraction and extends to higher tempera- 
tures. 

The role of the range of the attractive interactions in stabilizing fluid, self-assembled mem- 
branes qualitatively resembles the role it plays in stabilizing a fluid phase of simple molecules. 
If such molecules interact via a hard-core repulsion and a weak, but longer-ranged, attraction a 
fluid phase exists only if the range of the attraction is sufficiently large, roughly greater than 20% 
of the hard core diameter. Otherwise the fluid directly condenses into a solid [ I176| [T77| I178j . 

Another solvent-free model has been devised by Brannigan and co-workers [ 11791 [TSUI 11811 ITS^ 
I183j . The amphiphiles consist of five beads. The ffi'st bead (h) corresponds to the hydrophilic 
head, the second bead (i) is associated with the interface between hydrophilic and hydrophobic 
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Figure 2. Phase diagram resulting from Deserno's model as a function of rescaled temperature, 
ksT/e and range, Wc/cfzi, of the attraction between tail beads according to Eq. ((T^ . The areal 
density corresponds to zero tension. Each symbol corresponds to one simulation and identifies 
different bilayer phases: x - gel/crystalline, • - fluid, + - unstable. Lines are merely guides 
to the eye. The inset shows the pair-potential between tail lipids (solid line) and the purely 
repulsive head-head and head-tail interaction (dashed line). From Ref. [ I174j . 



entities, and the other three beads constitute the hydrophobic tail (t). The distance between 
neighboring beads along the amphiphile are constrained to a distance cr, which defines the length 
scale. A bending potential of the form Vbend(^) = ^bendcos^, with 5e < febend < lOe, tunes the 
geometrical shape of the molecules. A repulsive interaction 

12 

"^i-ep \ - \ '^i-ep = 0.4e, (15) 



Vrep(?' 



-rep 



is applied between all beads except intramolecular pairs separated by less than three bonds. 
Tail beads, and a tail and an interface bead, attract each other via a standard van-der-Waals 
attraction: 



att^ 



(16) 



Both, repulsion and attraction, are truncated at a distance, 2a. The interface beads, however, 
interact among each other via a longer-ranged potential 



Mnt(r) 



-Cint 



with 



Cint = 3e, 



(17) 



which is cut off at 3cr. This longer-ranged attraction, which acts only between the interface 
beads, is sufficient to stabilize fluid bilayer membranes at reduced temperature, ksT/e = 0.9. 

A different path to the development of generic coarse-grained models has been pursued by 
Groot [HHIIM], Smit [IIHElSLJIM. 189,290,191,, Shillcock [142,i92,i93,i94| and Mourit- 
sen [ I195| 11!??)] with their co-workers. These coarse-grained models utilize ultra-soft interactions 
in conjunction with a dissipative particle dynamics (DPD) thermostat [ 11971 [T^51 11991 1200l I201j . 
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Unlike the Langevin thermostat that adds random noise and friction to each particle, the DPD 
thermostat adds random noise and friction to each neighboring pair of particles. Thus, mo- 
mentum is locally conserved and hydrodynamic flow can be described. The use of ultra-soft 
potentials allows for rather large time steps for integrating the equation of motions (see Refs. [ 
121)21 I2U31 11961 12U4j for a detailed discussion). Their use can be justified by recognizing that 
the center of mass positions of the group of atoms that constitute a coarse-grained segment can 
overlap and their interaction is much softer than the harsh repulsions on the atomistic scale 
(cf. Sec. l2.2.3|l . This is a generic feature of coarse-grained models: the larger the length scale the 
weaker the interactions. By the same token, the density of the soft beads exceeds unity when 
measured in units of the particles' interaction radius. 

In DPD simulations, particles of type i and j (denoting water (w), head (h), glycerol-linking 
(e) and tail (t) ) interact via a very simplistic soft force of the form: 

F.,(r) = I - H fo-^^'^ (18) 

[ U for r > Tc 

where r is the distance vector between the particles' positions. The range of the interactions, 
Tc, between these soft beads sets the length scale. 

Originally, the DPD simulation scheme had been devised to simulate fluid flow, and a soft 
bead was thought of as a fluid volume comprising many molecules but still being macroscopically 
small. In the context of membrane simulations, a soft bead represents a rather small fluid volume 
that consists only of several molecular groups comprising the amphiphile. Often one identifies the 
range of interaction, Tc, with 1 nm, i.e., one tail bead corresponds to three or four methyl units. 
By the same token, a solvent bead corresponds to a small number of water molecules. (Attempts 
to map a single methyl unit onto a soft bead were rather unsuccessful [ 188 in reproducing the 
internal bilayer structure and resulted in a significant interdigitation of the apposed monolayers.) 
Typically the amphiphilic molecules consist of only a very small number of beads - one to three 
hydrophilic head beads and four to ten hydrophobic tail beads. The longer the amphiphiles the 
more stable and rigid the bilayer is. 

The strength of the interaction simultaneously describes the repulsion between unlike species, 
and the excluded volume of the coarse-grained beads which imparts a small compressibility to 
the liquid. The parameters of the model are tailored to reproduce key characteristics of the 
amphiphiles in solution (e.g., the compressibility of the solvent and the bilayer compressibility). 

A typical set of interaction strength is: 
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These values were originally proposed by Groot to parameterize the interactions of ionic surfac- 
tants [ .185) . Smit and co-workers [ I188j increased a^w from 15 to 25 to avoid very high densities 
in the bilayer's hydrophobic core. Different sets of parameters have been devised for double-tail 
lipids. Shillcock and Lipowsky [ 1142j used the set 
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Figure 3. Phase diagrams of model lipids 
as a function of head-head repulsion, ahh, 
and temperature, T* for various lipid ar- 
chitectures as indicated in the key. The dif- 
ferent bilayer phases are sketch in the bot- 
tom row: Lc ~ tilted subgel phase, Lj^ - 
gel phase, - tilted gel phase, - rip- 
ple phase. La - fluid phase. In the two 
lower phase diagrams "c" denotes a coex- 
istence region, of which the exact structure 
was difficult to determine. From Ref. [ I191j . 



Interactions along the amphiphile determine the molecular shape of the amphiphile. They 
include a harmonic bonding potential of the form [ I185j 

Vbond(?^) = ^bondr^, (21) 
or [UHS] 

^bond(r-) = /cbond(|r| -ro)2. (22) 
Additionally, bending potentials of the form [ Hi 



^bend — -^Kendif^ — ^o)^, (23) 



or [Eg 

Vbend = A:bend(l " COs(6l - 6*0)), (24) 

have been applied. Shillcock and Lipowsky [ I142j highlight the role of the bending potential on 
the internal structure of the bilayer membrane. The earlier DPD model by Groot [ I185j . and 
also self-consistent field models [ 1205] . describe the lipid tails as completely flexible which leads 
to a broad distribution of the tail ends throughout the bilayer. This interdigitation of the two 
apposing monolayers that form the membrane typically is not observed in biological membranes. 
A very large incompatibility between hydrophobic and hydrophilic entities would be required to 
stretch the monolayers sufficiently to reproduce the structure of a biological membrane. 

This example illustrates that bilayer structure and properties do sensitively depend on the 
model parameters of the minimal coarse-grained models. Since the potentials of these models 
are not directly related to molecular interactions, the parameters have to be adjusted so as to 
reproduce macroscopic observables. Although the minimal coarse-grained models are quite sim- 
ple, they have been shown to form several of the known bilayer phases: fluid, gel and crystalline 
[IMlIIEnilinnilMlIinHlIIIIlIIISl, as wen as more exotic ones, like the ripple phase [UHZHM]. 
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The different phases that can be obtained by the DPD model of Smit and co-workers are il- 
lustrated in Fig. ISI The phase diagram as well as the area compressibility modulus, bending 
stiffness and spontaneous curvature [ 11931 11831 12U9j have been explored as a function of model 
parameters. This ensures that one can choose parameters that make the coarse-grained model 
mimic, qualitatively, the behavior of an experimental system. 

2.2.3. Systematic coarse-graining: potential and limitations 

While the above minimal coarse-grained models are able to explore the generic features of 
the behavior of amphiphiles in solution, much recent interest has focused on "deriving" coarse- 
grained models for a specific chemical substance [ El I21()| ITU IT^ ITHl [T7| . The general feature 
of these systematic coarse-graining schemes is the attempt to utilize information about the 
microscopic structure, obtained by atomistic simulations, for instance, in order to parameterize 
the interactions between the coarse-grained entities. 

These techniques have originally been developed for polymer systems, and they have been 
extensively employed to describe quantitatively the structure and dynamics of polymer melts 
and solutions [ |H[ CHI CH EIH EEl ElSl EUl EEl ■ While the universal properties of amphiphilic 
systems only require a coarse-grained model to capture a few relevant interactions, the system- 
atic construction of such a model that quantitatively reproduces the behavior of the underlying 
microscopic system is a very ambitious task. From the onset, it is obvious that decimating 
the degrees of freedom will result in a loss of information, and will generate very complicated 
multi-body interactions [ I216j . The latter imparts additional complications on extracting ther- 
modynamic information from coarse-grained models. Commonly the multi-body interactions 
are, in turn, approximated by pairwise interactions in order to retain the computational effi- 
ciency of the coarse-grained representation, and these effective pairwise interactions depend on 
the thermodynamic state of the system (i.e., they depend on density or temperature, and they 
are different in different thermodynamic phases). 

The general principles of systematic coarse-graining consist in (i) choosing a set of key char- 
acteristics that the coarse-grained model shall reproduce, and (ii) determining the interactions 
between the coarse-grained degrees of freedom so as to reproduce these characteristics. The first 
step is the most crucial, and is guided by insight into the physics of the phenomena that one 
wants to investigate. Three qualitatively different properties can be distinguished - structural 
quantities, thermodynamic properties, and dynamic characteristics. 

Structural quantities are related to the geometry of the molecules on different length scales. In 
order to capture the specific details of the molecular architecture, a coarse-grained model should 
not only reproduce the overall size of the molecule (e.g., the end-to-end distance), but it should 
also include finer details of the molecular architecture (e.g., the stiffness along the molecular 
backbone, the bulkiness of the lipid's head, the location of double bonds in the hydrocarbon 
tails). It is essential to capture the rough features of the molecular geometry and its fluctuations, 
or ability to deform in response to its environment. 

In the ideal case the parameterization of a coarse-grained model is based on the detailed 
atomistic information about the molecular conformations under the same thermodynamic con- 
ditions, i.e., the same temperature and density. Then, one explicitly defines a mapping from a 
configuration of the atomistic model, {r}, onto a configuration of the coarse-grained model, {R}. 
Utilizing a large equilibrated sample of atomistic configurations, one can obtain the probability 
distribution of distances between coarse-grained entities. 

This procedure is illustrated in Fig. 0] for the case of a polybutadiene melt at 240 K and 
atmospheric pressure. The rich structure of local correlations is specific to the chosen system, 
but the qualitative features discussed in the following are born out by a wide class of systems 
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Figure 4. (a) Coarse-graining of the intramolecular segment-segment correlation function of 
polybutadiene (240 K). The numbers mark data with different degrees of the coarse-graining. 
For example, the coarse-graining degree n = 4 yields a chain molecule for which four united 
atoms are approximated by one effective segment. The bold curve is the correlation function 
calculated using the united-atom model (J-type peaks at r = 1.34 A, 1.5 A, 1.53 A arising from 
the CH=CH, CH2— CH and CH2— CH2 bonding distances are omitted), (b) Coarse-graining of 
the intermolecular segment-segment correlation function of polybutadiene at 240 K. For n < 8, 
the correlation function changes quantitatively. For large degree of the coarse-graining (e.g., 
n > 13), the correlation peaks and minima disappear and the correlation hole at r < 3A shrinks. 
The bold curve is the correlation function calculated using the united-atom model. From Ref. [ 
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and have been observed in a large variety of studies [ UHl [HI iEl HSl 11101 HHl HH UlSl 

ma nn ma usi mni na ma mni mni iisia mn ma msi ma msi- The two panels of the 

figure present the intramolecular and intermolecular pair correlation function as obtained from 
a simulation within a united atom model [ I236[ 123 7| I238[ I239j . The intramolecular correlations 
exhibit a rich structure on short length scales that mirrors both the correlations due to bond 
length and torsional interactions, as well as the delicate packing of the dense polymeric fluid. 
At large distances, those correlations have died away, and the intramolecular pair correlation 
function smoothly decays like fl'intral'^) ~ 1/^ as expected for a Gaussian polymer in a melt [ 
la . The intermolecular pair correlation function describes the probability of finding two coarse- 
grained segments on different molecules at a distance, r. It exhibits qualitatively a form that 
one also expects for a simple liquid. At small distances the correlation function vanishes - this 
distance characterizes the "thickness" of the polymer. There is a broad nearest-neighbor peak 
around 5 A, and there are a few further oscillations at larger distances. The long-range approach 
of the intermolecular correlation function to unity is dictated by the polymeric correlation hole 
effect and it is identical to the decay of the intramolecular correlations because the two cancel 
one another at length scales larger than the correlation length of density fluctuations due to the 
near incompressibility of the dense liquid [ ^ . 

One can use the explicit configurations of an atomistic simulation, and lump n successive 
carbon atoms along the backbone of the polymer into an effective coarse-grained segment. In 
this specific example, the location of the coarse-grained bead is taken to be the center of mass of 
its constituents, without accounting for the differences in molecular mass of CH and CH2 units. 
The correlations between the coarse-grained beads constructed in this manner for different levels, 
n, of coarse-graining are also depicted in Fig. |a Two general characteristics can be observed: 
First, the larger the degree of coarse-graining, the smoother are the intra- and intermolecular 
pair correlation functions. This behavior stems from the fact that the coarse-grained beads can 
partially overlap. This softening increases with the degree of coarse-graining, n. Second, the local 
structure is very sensitive to n, while the large distance behavior is not. In the specific example, 
constructing a coarse-grained bead from n = 3 segments yields a rather poor representation of 
the intramolecular and intermolecular correlations: The first peak of ^intra coincides with the 
minimum of the data for the united atom model and the first peak of Winter is too high. The 
value ra = 4 yields a better description. Increasing the value of n even further, one finds that 
the local structure gradually fades away and the beads become so soft that there is a significant 
probability that the coarse-grained beads overlap (c.f., inset of Fig. EJb). Thus there is an 
optimal degree of coarse-graining. Values of n which are too small lead to comparably harsh 
coarse-grained potentials which, in turn, give rise to packing effects that are not related to the 
structure of the underlying atomistic systems. Large values of n result in a loss of local structure 
and thereby of chemical specificity. 

The thermodynamic equilibrium properties of a coarse-grained model can be constructed 
formally via a partial trace, i.e., a summation over all microscopic configurations compatible 
with a fixed set of coarse-grained degrees of freedom. This procedure is in complete analogy with 
Renormalization Group calculations which have successfully been employed in polymer physics 
and critical phenomena [ 0101 EI- Let {r} denote the coordinates of the detailed microscopic 
model (e.g., the atom positions in a chemically realistic model) and let {R} denote the coarse- 
grained degrees of freedom. There is a mapping from the detailed degrees of freedom, {r}, onto 
the coarse-grained ones, R[{r}]. For instance, {R}, might be the center of mass of a small 
group of atoms or the location of a particular group of a lipid molecule. Let £'[{r}] denote the 
pairwise interactions of the microscopic system, so one can calculate an effective interaction, 
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C/[{R}] between the coarse-grained degrees of freedom via [ I216j 

e.J-Wl) = fvUr)]e.J-W)m-mr)]) (25) 



ksT J J \ ksT 

With this definition of the effective interaction, {/, the partition function of the microscopic 
system Z can be obtained according to 

and the probabihty of the coarse-grained degrees of freedom is identical to that generated by 
equihbrium configurations of the microscopic model. There are, however, two caveats: 

(i) U does not possess the typical characteristics of an interaction but rather those of a free 
energy. In particular, U depends on the thermodynamic state characterized by temperature, 
pressure, etc. Due account of this state dependence has to be taken if thermodynamic quantities 
are extracted. For example, the internal energy is calculated according to [ i240] : 



/P[{r}]S[{r}]exp 



E\\r}] 



jV[{r}] exp [-^^) 



^2 d In Z 



!V[{R}] {U[{R}] - T^U[{R}]) exp (- 



/P[{R}]exp(- ^ 



m{r}] 



(29) 



= (^U[{R}]-T-^U[{R}]j^^^^ (30) 

By the same token, the effective potential depends on density, p, and care has to be taken when 
calculating derivatives of thermodynamic quantities with respect to the number of particles or 
volume. For density-dependent central pair potentials, U{R,p,T, • • •) the pressure is given by [ 

dF 

= Pideal + Pvirial + ^ / d^i? ^Wl^IlllA g(^R^ p^T,---) (32) 



2 7 dp 

with Pidcs.l = ksTp and Pvirial = -y / d^RR g{R) (33) 

where the first terms are the familiar ideal gas term and the viral expression for the pressure 
which ignore the density dependence of the effective potential. The additional second term arises 
due to the dependence of the effective interaction, U, on the thermodynamic conditions under 
which it has been obtained [ 12411010] . g{R) is the pair correlation function of the coarse-grained 
fluid (cf. Eq. dSni)). 

(ii) Generally, ?7[{R}] cannot be decomposed into pairwise interactions but it consists of 
complicated many-body interactions. It is neither feasible to construct these many-body inter- 
actions numerically through the formal coarse-graining procedure outlined by Eq. ()25() (see Refs 
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[ 12421 12431 12441 1245j for explicit schemes for simple lattice models) , nor can they be utilized in 
an efficient computational model. 

In order to obtain a computationally tractable model, one has to approximate the effective, 
many-body interaction by pairwise potentials, and one often ignores the state dependence of the 
effective interaction. To this end, one seeks an approximation of the form 

[/[{R}, p,T,...]^J2J2 ^(^^ - ^j) (34) 

i j<i 

The detailed form of V is unknown. Often one utilizes a generic form that contains several 
free parameters. These parameters are adjusted to the specific system (see below). Moreover, 
one requires that even after ignoring the state dependence of the interactions, the model still 
reproduces thermodynamic properties such as the equation of state, the interfacial tension, etc. 
It is this uncontrolled approximation in the construction of coarse-grained models that requires 
physical insight. 

Once the general type of the interactions to be included and the key characteristics the coarse- 
grained model should reproduce have been chosen, the second step of constructing the detailed 
form of the coarse-grained interactions, V, is conceptually straightforward but tedious. As a 
first approximation one utilizes the potential of mean force, V^f obtained via the pair correlation 
function of coarse-grained degrees of freedom: 

5(R.,)^exp(-^^^|^) (35) 
with Rjj = Rj — Rj , and where the pair correlation function is given by 

g{R) =V{6iK- {R,[{r}] " RjlM]))) (36) 

through the microscopic configurations. 

By definition this construction procedure works well if the coarse-grained degrees of freedom 
are dilute, but fails if higher order correlations become important (condensed phase effect). Nev- 
ertheless this potential of mean force serves as an initial guess for more sophisticated procedures. 
One strategy consists in using the Ornstein-Zernike equation 

/i(Ri2) = c(Ri2) +pJdK3 /i(R23)c(Ri3) (37) 

to obtain the direct correlation function, c(R) from the total correlation function, /i(R) = 
g(R-) — 1. Then, one can use an approximate closure relation to express the direct correlation 
function in terms of the interparticle potential. Liquid state theory provides guidance in choosing 
accurate closures. An example is the hypernetted chain closure (HNC) [ 1232j 

^ = -ln5(R)+MR)-c(R) (38) 

which differs from the potential of mean force in Eq. (|35() by the last two terms. 

This scheme provides a systematic, yet approximate, approach for constructing the effec- 
tive pair interactions. In polymer physics, one often employs a numerical method - iterative 
Boltzmann inversion [ 12461 12471 12311 123U1 1232j . One starts with an initial guess for the pair 
interactions, Vq. This pair interaction is used in a computer simulation of the coarse-grained 
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system, and the correlation function, go, is monitored. In the next step of the iteration, i = 1, 
the potential is improved according to 



where g{Ti) is the reference pair correlation function obtained from the microscopic system. 
The procedure is continued until the coarse-grained simulations with the potential Vi reproduce 
the reference pair correlation function with sufficient accuracy. Other iterative scheme, e.g., 
self-adjusting Wang-Landau-type methods [ I248j . can be envisioned as well. 

In addition to these structural quantities the potential parameters of the coarse-grained repre- 
sentation are tailored to reproduce thermodynamic properties such as pressure, interface tension 
or elastic moduli. Some quantities (like pressure or membrane tension) can be included in the 
iterative Boltzmann inversion scheme utilizing constant pressure or tension simulations [ l23Uj . 
Alternatively, non-bonded interactions that do not strongly impact the molecular structure are 
adjusted to reproduce the desired thermodynamic properties. 

While the mapping of equilibrium properties between atomistic and coarse-grained models is 
conceptually well understood, the mapping of dynamical properties might prove an even bigger 
challenge. The softer interactions in coarse-grained models typically speed up the dynamics on 
large length scales. This is a major computational benefit of coarse-graining. However the speed 
up is not uniform, but rather depends upon the specific dynamic property. For instance, the ratio 
of time constants between lateral diffusion of the lipids' center of mass, the structural relaxation 
of the lipid conformations and the rotational diffusion around the bilayer's normal differ between 
the coarse-grained model and atomistic simulations [E]- Several factors contribute to the non- 
uniform speed-up: 

The mobility of lipid head groups is coupled to the structural rearrangement of water in 
their vicinity, but the local structure (e.g., hydrogen bonding) is not faithfully reproduced by 
coarse-grained models. Likewise, the elimination of degrees of freedom of the hydrophobic 
tails (C-C bond stretching, bond angle and torsional potentials) in the course of the coarse- 
graining procedure speeds up the structural relaxation of the lipid tails. Moreover, the neglect of 
topological constraints in the molecules' dynamics (non-crossability) by very soft potentials (e.g., 
DPD models) will selectively speed up the dynamics of extended molecules compared to e.g., 
small solvent molecules. These effects have quite different origins and, in general, will result in 
different speed-up factors. Moreover, electrostatic interactions omitted in coarse-grained models 
might be not essential for structural properties, but still may significantly alter the dynamics 
[ I249j . Likewise, integrating out the solvent molecules results in strong attractions between 
like units of the amphiphilic molecules. This might dramatically influence the dynamics (see 
e.g., Ref. [ I250j for a careful discussion of the role of explicit solvent on the dynamics of the 
collapse of a polymer chain in a bad solvent). Therefore, one cannot expect that a single 
scale factor relates the characteristic times of different motions and dynamic processes between 
atomistic and coarse-grained models [ I13| I15| I251j . Thus particular care has to be exerted when 
non-equilibrium process are to be described by coarse-grained models. Nevertheless, a natural 
and common identification of the time scale in coarse-grained models consists in matching the 
characteristic time of the slowest process of relevance, e.g., the lateral diffusion of a amphiphilic 
molecule in the bilayer. 

2.3. Coarse-grained field-theoretic models and molecular field theories 

An alternative approach to a particle-based description has its roots in statistical field theory 
[ I252| I253j . In this framework one operates with collective variables, such as density and in- 
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teraction fields instead of tlie coordinates of the individual molecules. The formal advantage of 
this approach is that it allows one to decouple many-body (particle-particle) interactions and to 
replace them with one-body (particle-field) interactions. This "simplification" comes at a price 
- to calculate the partition function one has to perform thermodynamic averages over the field 
variables. Nevertheless, under suitable conditions, one can introduce controlled approximations 
that make the field-based approaches tractable. 

Formally, the particle-based canonical partition function of a system of n particles, occupying 
volume V and interacting through a (multi-body) potential U{{r}) 

n „ 

Z oc n / d=^r,- e-^(ir})/fcsT^ (40) 
can be transformed into the following field-theoretic partition function: 

Z oc / Vwo, f V<Pp e-^['"-<^^l/'=^^, (41) 



where '^[wa, (^/3] is usually called an effective action or Hamiltonian of the system, and it depends 
on density and interaction field variables, (/>/3(r) and Wa{r). Here, the indices a and (3 run over 
different types of microscopic units (monomers, segments), and the functional integration JP/ 
is performed over all spatial realizations of the field variables. 

Despite the fact that the field-theoretic transformation replaces one very hard problem, 
Eq. H40() . with another, Eq. H41|). in this representation it is simple to implement the mean- 
field approximation which, in polymer systems, is quite accurate [ 12541 I255j . Physically, the 
mean-field approximation amounts to the assumption that the statistics of a system is domi- 
nated by a particular field configuration, so that the partition function is also dominated by a 
single term and takes the simple, approximate form 

Z oc e-«[<'<^51/'=s^, (42) 

where the dominant field configurations, and (j)*^, are the stationary (and, actually, saddle 
point) configurations of the effective Hamiltonian: 



6n 



= (43) 

This approximation is usually not very accurate for systems composed of low-molecular weight 
species with short-ranged interactions, but it becomes more so for particles with long-ranged 
interactions and for high-molecular weight polymers at high concentrations. Depending on 
the system of interest, one can frequently quantify deviations away from the mean-field solution 
caused by long- wavelength fiuctuations, and define the corresponding Ginzburg parameter [ [2541 
[255 . In polymer systems with short-ranged segment interactions, it is usually the dimensionless 
overlap parameter, or invariant degree of polymerization, J\f = (pRl/N)'^, which characterizes 
the number of polymer chains interacting directly with a given chain. Here Re is the molecules' 
end-to-end distance, N denotes the number of monomeric units (segments) per molecule and, p 
is the segment number density. This parameter can be compared to that in the case of simple 
lattice models. In the latter systems, the mean-field approximation usually becomes accurate at 
large enough site coordination number, commonly achieved either at high spatial dimensionality 
or for long-ranged interaction potential. In polymer systems, however, high effective coordination 
numbers can be achieved by increasing the polymerization index at a constant segment density. 

There are two basic issues that have to be addressed within any mean-field approach: (i) 
calculation of the average molecular-field acting on a given molecule due to a density distribution 



Coarse-grained models for biological and synthetic membranes 



27 



of all other molecules present in the system, (ii) evaluation of a single-molecule partition function 
in the presence of this field, which can be formally written as: 



Here w{s) is the mean-molecular field felt by the chain in microscopic state s. These states 
have o priori statistical weights gs , that are solely determined by intra-chain degrees of freedom 
and characterize the molecular architecture in a spatially homogeneous system. Ultimately, the 
two ingredients of a mean-field approach have to be solved self-consistently for all the species, 
therefore the name self- consistent field (SCF) theory. The variety of the SCF approaches in the 
literature can be classified by the types of approximations used for solving the two problems 
above. Below, we describe a number of such approaches as applied to modeling of bilayer 
membranes. 

2.3.1. Anchored chain models 

One of the first molecular- field theories of a lipid bilayer was proposed by Marcelja [ 12561 1257] . 
In this model, lipid head groups and solvent molecules are modeled by a simple boundary to 
which lipid tails are anchored at a given areal density. Intra-chain degrees of freedom are 
sampled using Flory's Rotational Isomeric State (RIS) model [ 12581 I259j and the segments 
interact through an anisotropic aligning potential of Maier-Saupe kind: 



where 0i2 is the angle between segment axes. This potential is borrowed from the theory 
of nematic liquid crystals, and in the context of lipid packing the aligning tendency between 
segments stems from the fact that all- trans chains pack at a higher density and, hence, experience 
stronger dispersive van der Waals attractions. Further, being better aligned, they are subject to 
fewer hard core repulsions. In the mean-molecular field approximation the two-body interactions 
in Eq. ()45() are replaced by a one-body interaction with the scalar order parameter field: 



where denotes now the angle between the segment and the average axis of orientation. This 
relatively simple theory is capable of a quantitative (albeit phenomenological) description of 
experimental data on the liquid-gel ordering transition thermodynamics. 

Gruen [ 12601 12611 12621 I263j extended Marcelja's model by replacing the phenomenological 
Maier-Saupe ordering interaction with a local incompressibility condition in the lipid tail region 
imposed through an inhomogeneous Lagrange multiplier field, 'k{z), coupled to the local micro- 
scopic segmental density (j){z) = Ylf=i^{^ ~ ^i)) where z is the coordinate across the laterally 
homogeneous bilayer and the sum runs over all the positions Zi of the N units per molecule. 
Results derived with this model show semi-quantitative agreement with results of molecular 
simulations both for thermodynamics as well as for the structural properties of bilayers. Minor 
modifications of this approach allow the study of properties of cylindrical and spherical micelles 

Similar ideas of imposing packing constraints on anchored chain systems have been explored 
by a number of researchers. Differences between these methods are mostly technical. Dill et al. [ 
[267, 268_ and Scheutjens and co-workers [ I269j use a lattice to describe both chain conformations 
as well as molecular field configurations, which can be treated with efficient transfer matrix 
methods originally proposed by Kramers and Wannier for solving the Ising model. Ben-Shaul 
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et al. [ I27UI 12711 12721 1273j enumerate a very large set of the most probable chain configurations 
off-lattice and use iterative techniques to reach self-consistency. For a more detailed exposition 
of these ideas we refer the reader to excellent reviews of these works [ 12681 127Uj . 

Recently EUiott et al. [EZHEZSl has extended the mean-field approach to study mixtures of 
saturated and unsaturated lipids and cholesterol. It has been shown that the incompressibility 
(or packing) constraint is not sufficient to provide a quantitative description of both liquid and gel 
lipid membranes and the corresponding main-chain transition. This deficiency can be eliminated 
by including additional explicit orientational interaction, similar in spirit to the original ideas of 
Marcelja and Gruen. The following form of pairwise interaction potential has been proposed: 

V{u, u') ^ V{u • c, u' • c) « - J/(u • c)/(u' • c), (47) 

where u and u' are the local orientations of CH2 groups, c is the bilayer normal, and J charac- 
terizes the coupling strength. In the absence of head group orientational interactions, the chain 
tails on average will align parallel to the bilayer normal, which implies that f{x) should be a 
monotonically decreasing function of its argument. Elliott et al. [ I274j utilize the generic form 

9777 -I- 1 

/(cose) = — ^(cos^e)"^ « mexp(-me2), (48) 

where is the angle between the local chain orientation u and the bilayer normal c. The 
parameter m controls the width of the orientational interaction and the case m = 1 is comparable 
to the original Maier-Saupe interaction function used by Marcelja [ 1276j . Elliot et al. [,274, take 
m = 18. The higher power dependence was first proposed by Gruen [ 1262| 1263] . 

The experimental value of the main-chain transition temperature of a single component DPPC 
system is used to fit the strength, J, of the orientational interaction whereas the power, m, is 
estimated from the average alignment of an acyl tail in the bilayer. The same set of param- 
eters can be used to describe membranes containing unsaturated dioleoylphosphatidylcholine 
(DOPC) lipids. This theory successfully predicts the absence of the liquid-gel transition in 
DOPC membranes due to disordering effect of the unsaturated double bond. Binary dipalmi- 
toyl phosphatidylcholine (DPPC)/DOPC model bilayers exhibit phase coexistence between a 
saturated-lipid-rich gel phase and an unsaturated-lipid-rich liquid crystalline phase with the 
transition temperature below that of the pure saturated-lipid system. An extended model of 
a ternary system containing saturated and unsaturated lipids and cholesterol has been used to 
provide strong evidence for the existence of a liquid-ordered to liquid-disordered phase transition 

Despite the fact that the anchored chain models provide a detailed description of chain- 
packing in both gel and liquid lipid bilayers, their description of the interfacial polar-head region 
is very rudimentary. It is either completely ignored in favor of fixing the lipid areal density or 
treated essentially phenomenologically. The inequivalence of tail, head and solvent segments 
effectively makes the bilayers pre- assembled structures, as opposed to self-assembled ones. Since 
this limitation does not come from the mean-field approximation, it can be completely eliminated 
within a more general self-consistent field-theoretic framework. 

2.3.2. Self-assembled membrane models 

The self-consistent theory of polymer systems, both for solutions and melts of long-chain 
molecules, was formulated by Edwards [ 1277j . Helfand and Tagami [ 278, 279, 280^ pioneered 
the use of this approach for the study of multicomponent polymer systems with diffuse interfaces. 
Since then the self-consistent field theory has been utilized for a wide range of polymeric systems 
[EOni ESD M EHSl EHi ESS ESS M, notably it has been highly successful in studying self- 
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assembly of block copolymers [ 12861 0571 1288j and copolymer-homopolymer mixtures [ 128911^^01 

Leermakers et al. used similar ideas to study self-assembled bilayers composed of polymer 
chains with random- walk and RIS (Rotational Isomeric State) statistics [ |205] . For the sake 
of computational convenience, chain configurations are restricted to a tetragonal lattice aligned 
with the bilayer interface. The segment-segment dispersive interactions are modeled by a set 
of Flory-Huggins parameters: XAW for tail segment /water, xbw ^or head group/water and 
Xab{^ Xaw) tail segment/head group interactions. The harshly repulsive interactions at 
short distances are taken into account by imposing the common incompressibility condition 
(j)Aiz) + (j)B{z) + 4>w{z) = 1, where (t>A,B,w{z) are local volume fractions of the tail, head, and 
water segments, respectively. The authors use two types of chain statistics to describe the in- 
tramolecular conformations. The first one is a random walk on the tetragonal lattice, which 
allows neighboring bonds to overlap and does not distinguish between the trans and gauche con- 
formations. The second more elaborate model, based on the RIS scheme, takes into account the 
small energy difference between the local bond configurations and also forbids direct backfolding 
of the chain onto itself. In both approaches the chain partition function and density distribu- 
tions can be evaluated using straightforward matrix algebra. Nevertheless, for computational 
tractability only laterally homogeneous structures are considered. 

This fully self-consistent framework is capable of describing stable, tensionless, self-assembled 
bilayers, provided the interaction parameters are within certain limits. The random-chain and 
the RIS-chain models both result in membranes with qualitatively similar segment distribu- 
tions and thermodynamic properties. Quantitatively, however, this approach underestimates the 
experimentally measured membrane thickness by about 50%, suffers from an excessive chain- 
disorder and fails to predict the gel-liquid main-chain transition. Despite these deficiencies, this 
simple framework is very fiexible and efficient and it has been extended to studies of bilayer 
vesicles, membrane elasticity, bilayer-bilayer interactions and effects of surfactants on bilayer 
structure and stability [ I294j . 

The absence of the main-chain transition in this model has been traced to the inability of 
the simple incompressibility constraint to properly describe anisotropic packing of the chains. 
To cure this problem, Leermakers et al. formulated an anisotropic SCF theory [ I295j . The 
idea is based on a generalization of Flory's and Di Marzio's treatment of rigid rod statistics 
in anisotropic states and results in an effective entropic aligning interaction between chains 
(cf. Elliott et al. [ I274j ). The modified theory predicts a first order main-chain transition and 
the existence of two gel phases: one with intercalated monolayers in the dilute regime, and the 
other with almost non-overlapping monolayers in the concentrated regime. 

With straightforward modifications, this theory was used to investigate mechanical properties 
of the lecithin bilayer membranes. To calculate the bending moduli of the lecithin vesicles, 
Leermakers et al. [ I294| proposed to use spherically curved lattices that can accommodate either 
lattice random- walks or RIS chains. The fact that the underlying lattice is "commensurate" 
with the segment structure of the chains, leads to strong lattice artifacts. These numerical 
artifacts give rise to strong variations of both structural and thermodynamic properties as a 
function of vesicle radius. To resolve this problem, the authors propose to consider only vesicles 
with radii that are commensurate with the underlying lattice. Despite the fact that the model 
then provides a reasonable account of the stability and structure of the curved, pure lecithin 
bilayers, this approach has been criticized in a later work (see below). Bending moduli of 
multicomponent vesicles are predicted to be smaller than those of pure vesicles composed of 
either of the components. This prediction is in agreement with the earlier studies on anchored 
chains grafted to curved interfaces [ I296j . 
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Oversteegen et al. [ I297j studied in detail the dependence of elastic constants on microscopic 
interaction parameters. In particular, they found that the bending modulus, Kb, strongly depends 
on the strength of interactions between the head groups and the solvent. As the strength of head 
group hydration decreases, so does the bending modulus, apparently vanishing for sufficiently 
weak head group/solvent interaction parameters. Interestingly, under the same conditions the 
saddle splay modulus, R, progressively becomes less negative and eventually adopts positive 
values. This favors formation of saddle- like structures and should lead to the instability of a flat 
bilayer. The authors also carefully investigate the spherical lattice artifacts mentioned above 
and conclude that there is actually no completely consistent way to eliminate them, despite 
earlier claims to the contrary. Nevertheless, the authors argue that the lattice model may still 
prove to be a valuable tool to to make qualitative predictions. The same approach was used 
[ I298| to study elastic properties of surfactant monolayers. It was found that the monolayer 
spontaneous curvature does not vanish, as in the case of the bilayers, and strongly depends on 
molecular architecture. The monolayer bending modulus was found to be essentially half of that 
of the corresponding bilayer, in agreement with previous studies which used an anchored chain 
model [Uni- 

Meijer et al. [ ISOOj used extensions of the lattice mean-field model to study inclusions of 
foreign molecules in bilayer membranes utilizing a detailed description of both rigid and flexible 
molecules with, and without, electrostatic interactions. To handle the rigid inclusions whose 
segment structure is not necessarily commensurate with the underlying tetragonal lattice (used to 
generate RIS hydrocarbon chain configurations), the authors must make several approximations 
to preserve efficiency of the propagator scheme. This approach is implemented in a computer 
program called GOLIATH. Partitioning coefficients for a wide range of molecules into bilayer 
membranes were predicted in reasonable agreement with experimental data. Detailed theoretical 
information on the density and orientational distributions of these inclusions allowed the authors 
to conclude that it is possible to design isomeric molecules that can have roughly the same 
partitioning coefficient, but would differ greatly with respect to the position and orientation 
of the inclusion in the membrane. This could be potentially of use in designing membrane- 
incorporated drug systems. 

In a more recent work, Leermakers et al. refined their SCF approach by including several new 
features [ l3Ulj . First, the introduction of different segment volumes for CII2 and CH3 groups 
counteracted the interdigitation of hydrocarbon chains from apposing monolayers. Second, by 
allowing the formation of clusters of water molecules [ 1302 j , they were able to prevent excessive 
penetration of the water molecules into the membrane's hydrophobic region. Third, polariz- 
ability of chain segments was included for thermodynamic consistency. The predictions of this 
detailed approach compare favorably with all-atom MD simulations by the same authors [ I303j 
of SOPC and SDPC bilayers. This validates the SCF parameter set, which can then be used to 
investigate other membrane properties. 

Whitmore et al. extended the lattice, self-consistent, anisotropic field theory to compressible, 
fully hydrated, phospholipid membranes [ l3U4j . The authors assumed that the lipid headgroups 
are strongly confined both positionally and orientationally throughout the temperature range 
of interest. This approximation considerably simplifies the numerical solution of the set of self- 
consistent equations, and emphasizes the bond density and orientation in the hydrocarbon region 
of the bilayer. The assumption, however, is in disagreement with the the rather wide interfaces 
observed by Leermakers et al. in a similar model bilayer system [ I301j . Similarly, the solvent 
molecules (water) were excluded from the hydrocarbon layer. To determine the microscopic 
interaction parameters the authors used the main-chain transition temperature for DPPE, its 
dependence on the chain length, and the change in thickness of the hydrocarbon layer at the 
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transition. The so-parameterized model exhibits a first order phase transition between the highly 
ordered gel phase, L/j, to the partially disordered liquid crystalline phase, La- It also provides 
a reasonable account of equilibrium changes in the bilayer thickness and hydrocarbon density 
across the main chain transition. The predicted increase of the transition temperature with the 
hydrocarbon chain length is in qualitative agreement with experimental observation, but is still 
weaker by a factor of 3. The authors attribute this inaccuracy to the neglect of orientational 
effects of water on the gel-to-liquid transition. 

In an alternative approach, Miiller and Schick [ 1305] used an off-lattice representation of the 
field theory, and obtained the single-chain partition function via a partial enumeration [ l3U6j over 
a large set of molecular conformations of a lipid chain with the RIS statistics. The model takes 
into account an explicit solvent, which permits the study of thermotropic, as well as lyotropic, 
phase transitions between lamellar, L^, and inverted hexagonal, Hu, phases. In agreement with 
experiment, it was found that the transition from Hu to occurs upon increasing the solvent 
content or decreasing the temperature. It was also confirmed that an increase in the length 
of the hydrocarbon tail or decrease in the head-group volume stabilizes the inverted hexagonal 
phase. 

In a series of papers, Schick and coworkers used a flexible chain model to study bilayer mem- 
branes. This model is known to be incapable of producing the main-chain transition to the gel 
phase, Ljs, due to the excessive entropy of a flexible chain. Nevertheless, in the liquid crystalline 
phase. La, it is expected to account reasonably well for the competition between hydrophobic 
forces and chain stretching. In particular, Netz and Schick [ l3U7j used a standard flexible diblock 
copolymer model in the presence of a solvent, modeled by a flexible hydrophilic homopolymer, 
to study structure and stability of the fluid bilayer phase. It was found that under external com- 
pressive stress the bilayer thickness decreases until the bilayer is in thermodynamic coexistence 
with either the disordered or catenoid lamellar (hexagonally perforated lamellar) phases. This 
result was used to describe membrane rupture. It was found that the relative areal expansion 
at the onset of rupture is about 2%, which is in agreement with experiments on lipid systems. 
Pores, which appear in the hexagonally-perforated lamellar phase, formed reversibly with stress 
and were slightly hydrophobic. 

Li and Schick [ 1308] used a more detailed microscopic model of a lipid with a charged head- 
group, two flexible hydrophobic tails, and a neutral solvent with counter ions. The hydrophobic 
interactions are described through the usual contact interactions. First, the model is applied 
to a neutral, anhydrous, lipid and is shown to describe reasonably well the transitions between 
lamellar, hexagonal and cubic phases. With the addition of a solvent, the system exhibits a 
re-entrant Hu La ^ Hu transition, a peculiar feature observed experimentally in the di- 
oleoylphosphatidylethanolamine (DOPE)-water system. Addition of a negative charge on the 
headgroup leads to an effective increase of the head-group volume and drives the transition from 
the Hu to the La phase, also in agreement with experiments. In a second paper, Li and Schick 
[ I3f)9j extend their approach to study mixtures of cationic, anionic and neutral lipids. This 
model provides a detailed description of the pif-tunable La — Hu phase transition in mixed lipid 
systems, which could play an important role in liposomal drug delivery systems [ I31flj . 

We used the SCF approach with the Gaussian chain model to elucidate the evolution and free 
energy costs of intermediates involved in the fusion of bilayer membranes [ l311ll^T^ . This study 
goes beyond the previous work on laterally homogeneous bilayers and provides a description of 
complicated structural rearrangements that occur during the fusion reaction. Detailed results 
of this approach are presented in Sec 13.41 Here we mention important modiflcations necessary 
to study such complicated problems within the field-theoretic framework. 

In practice, the SCF theory in its common implementation is only capable of identifying ther- 



Coarse-grained models for biological and synthetic membranes 



32 



modynamic, locally stable configurations of a system. Such structures can correspond to stable 
configurations of a system or metastable intermediate states encountered along the transforma- 
tion of the system towards a stable structure. In studying activated processes (e.g., rupture, 
fusion or fission of membranes), however, one should also be able to describe properties of transi- 
tion states which correspond to saddle points of the SCF free energy functional. Unfortunately, 
finding a saddle point of a functional poses a serious numerical problem. To address this issue 
we propose the use of reaction coordinate constraints [ I293j . In some cases, using either physical 
intuition, qualitative input from numerical simulations, or SCF fiuctuation calculations [ I287j . 
one can identify the unstable directions of a saddle point and stabilize them by applying a suit- 
able constraint. Such a constraint, for instance, can be a restriction on the symmetry of the 
segment density distribution, or a direct constraint on the field variables. 

As an example, consider nucleation of a hole in a bilayer membrane. First, we confine the 
possible solutions to be axially symmetric, and to have a mirror symmetry with respect to the 
bilayer mid-plane. We further constrain the possible segmental density field configurations by 
requiring that in the bilayer mid-plane, the interface between the hydrophilic and hydrophobic 
segments be localized along a circle of some prescribed radius, which plays the role of a reaction 
coordinate. This constraint can be imposed via a Lagrange multiplier field and results in a 
minimal modification of the usual set of SCF equations. In general, for an arbitrary value of the 
constraint, the value of the corresponding Lagrange multiplier field is non-zero, which means 
that the corresponding configuration of the system is not a solution of the original unconstrained 
problem. Therefore the details as to how the constraint is applied matters. Nevertheless, exactly 
at the points where the Lagrange multiplier field vanishes, the solutions coincide with those of 
the unconstrained system, and it is an extremum, a minimum, of the free energy along the 
reaction coordinate path, i.e., these points are intermediate states. Furthermore, an extremum 
which is a maximum of the free energy along the reaction coordinate path, corresponds to a 
saddle point of the unconstrained problem, i.e. to the transition states of the system along the 
considered reaction coordinate. This approach provides a wealth of additional information that 
cannot be accessed by the standard, non-constrained, calculations. 

Yet another approach that uses the collective density degrees of freedom as the basic degrees 
of freedom is the Density Functional Theory (DFT) . This formalism has its roots in the quantum 
many-body theory, where it is widely used for calculating electronic structure of many-atomic 
systems. Similar ideas have been applied to the description of simple classical fluids, in particular 
to the studies of vapor-liquid and liquid-solid interfaces, nucleation of new phases (see Ref. [ I313j 
for an application to hole nucleation in amphiphilic bilayers), and solvation of extended solutes. 
One of the major advantages of this approach is the systematic incorporation of liquid state 
packing effects on the level of segments. Chandler and co-workers [ I314j have developed an 
approach to study properties of inhomogeneous molecular systems, which combines molecular 
simulations and DFT. Despite the fact that the DFT is formally an exact framework it is limited 
by the accuracy of the approximate density functional. The common approach to construct an 
approximate density functional relies on mean- molecular field ideas, similar to those utilized in 
the SCF framework [ I255j . In particular, the interacting many-particle problem is reduced to 
that of a reference system consisting of non-interacting particles in an effective external field 
chosen so that it results in distributions equivalent to those in the original interacting system. 
Not surprisingly, in complete analogy with the SCF theory, the problem can be divided into 
two coupled sub-problems. First, one determines the mean- molecular fields produced by a given 
density distributions. Second, one evaluates the thermal averages of the density operators in the 
presence of these fields. The major difference between the DFT and SCF theory approaches lies 
in the microscopic details utilized to obtain the expressions for the mean-molecular fields. In 
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particular, in the DFT framework, the fields are customarily obtained within a linear response 
theory of a bulk uniform fluid to a weak perturbation. For example [ I315| I316j . the mean- 
molecular field Wa (r) acting on a segment of type a can be approximated by 

w.ir) = y„(r) -J2f dr'c^f (r - r') [p^ir') - p}"'^] , (49) 

where ppir') and p^^^^ are the local density and the bulk density of site type /?, respectively, Va{v) 
denotes the "bare" external fields (e.g., interactions with confining surfaces) and c^^^'^(r — r') 
is the bulk direct correlation function between segments of the types a and (3. The direct 
correlation function provides information about microscopic pair correlations in the system. 
In the mesoscopic SCF approaches discussed above, these microscopic correlations are usually 
treated on more phenomenological grounds and replaced by a simple contact interaction. The 
additional microscopic detail of the DFT approaches comes with the additional cost of calculating 
the direct correlation functions of a molecular fluid. This can be obtained either from simulations 
(see Eq. I37|) or from various integral equation approaches. Despite the much more elaborate 
microscopic input to the DFT framework, it is a priori not clear whether the density functional 
itself, obtained in the linear response approximation, is accurate in the strongly inhomogeneous 
systems such as self-assembled bilayer membranes. More sophisticated relations between flelds 
and densities have been utilized, e.g., to capture the equation of state of polymers in bad solvents 
[EmEIHlEin] or the local structure of the segment fluid [ElBES]- 

Frink et al. applied the DFT approach to study coarse-grained lipid bilayers [ I322j using 
a freely jointed chain model for the lipid tails. The direct correlation functions of the bulk 
disordered system were obtained from the polymer reference interaction site model (P-RISM) 
[ I323| I324| I325| I326j in conjunction with the Percus-Yevick closure. Additional single chain 
Monte Carlo simulations had to be performed to obtain the intramolecular structure factor, 
which is a necessary input to the P-RISM approach. The bilayer solutions obtained with this 
DFT framework have structural properties that are in a reasonable agreement with experimental 
DPPC data. Interestingly, the model even predicts a highly ordered bilayer phase at a sufficiently 
low temperature, which has an apparent resemblance to the ordered gel phase, Lp. Nevertheless, 
the disorder-order transition appears to be continuous, in contradiction with the experimentally 
observed first-order main chain transition. The reason for this discrepancy presumably lies 
in the high entropy of the fiexible chain model. A more accurate treatment would require 
the inclusion of chain stiffness to properly describe the main chain transition. Furthermore, 
in their comparison of this DFT approach to the MD simulations of the same model system, 
Frischknecht et al. [ I327j found no ordering transition in the latter at all. Therefore, despite the 
overall agreement with the MD simulations at relatively high temperatures, it seems that the 
DFT approach can become qualitatively inaccurate at low temperatures. 

2.4. Phenomenological Hamiltonians 

While the coarse-grained models discussed above retain the notion of lipid molecules, phe- 
nomenological Hamiltonians take the coarse-graining idea one step further and describe the 
configuration by the collective density of hydrophilic and hydrophobic entities or the spatial po- 
sition of interfaces. These approaches address the universal properties of amphiphilic systems, 
and are stripped down to essential ingredients. This permits a study of wide parameter ranges, 
and a systematic exploration of universal properties and mechanisms. The form of the Hamilto- 
nian is usually dictated by symmetry considerations. The few parameters of the Hamiltonian can 
be established through comparison of simple properties which result from the calculation, such 
as interface tension or fiuctuation spectra, either with experiments or with simulation results 
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of more detailed models. Two major approaches can be distinguished: (i) Ginzburg-Landau 
type models which describe the system via collective densities of hydrophilic and hydrophobic 
entities, (ii) Helfrich Hamiltonians which conceive of the system as an assembly of infinitely thin 
interfaces characterized by their elastic constants. We briefly discuss both approaches in turn. 

2.4.1. Ginzburg-Landau type models 

Ginzburg-Landau type models resemble the field-theoretic approach summarized above be- 
cause they describe the system via collective densities. Their staring point is a free energy 
functional that, in the simplest case, describes the free energy of a system with a spatially 
varying distribution of hydrophilic particles. Within the Random Phase Approximation [!^ one 
can derive such a Hamiltonian from the field-theoretic formulation of a microscopic model. This 
technique is useful for rationalizing the general form of the Hamiltonian, but the assumption of 
weak segregation between hydrophilic and hydrophobic entities inherent in the Random Phase 
Approximation is not valid in biological systems. Hence the approach is often quantitatively 
unreliable. Typically, the Hamiltonian is constructed via an expansion in terms of an order 
parameter such as the difference between hydrophilic and hydrophobic densities, (j){r), and its 
spatial derivatives, V0. Only the lowest terms consistent with isotropy of space are retained; 



The function determines the phases of the homogeneous system. Generically, it exhibits 
two or more minima. The function / can take the form of a polynomial, as in the original 
Ginzburg-Landau model of phase separation, or the form of a piecewise parabolic function of 4> 
[ I328j . The second term determines the free energy costs of interfaces. In phases characterized 
by an extensive amount of interface, such as a microemulsion, the coefficient in this term must 
become negative, at least for some composition range, in order to generate these stable interfaces. 
The third term with c > guarantees that spatial variation in the order parameter will not grow 
indefinitely, and therefore ensures the stability of the system. 

Although simple, this form of Hamiltonian suffices to reproduce many structures observed in 
lipid-water mixtures including disordered structures like microemulsions and complex periodic 
phases like the gyroid morphology [ I329j . Typically, the spatial derivatives are discretized on a 
regular cubic lattice. Then, the densities at the different lattice nodes are the degrees of freedom 
and the derivatives give rise to nearest and next-nearest neighbor interactions. The statistical 
mechanics of the model can be studied with Monte Carlo simulations [ l33Uj which includes 
thermal fluctuations. The fluctuations are important, e.g., for the formation of microemulsions. 
Alternatively, one seeks spatially modulated structures that minimize the Hamiltonian, and 
thereby one ignores fluctuations [ I329j . 

The attractive features of this type of model is the computational ease with which fluctuations 
or complex spatial structures can be investigated. Rather large systems can be studied over a 
wide parameter range, and often one can obtain analytic solutions in various asymptotic regimes. 
The connection between the model parameters and biological systems often remains qualitative, 
and we shall not discuss this approach further but rather refer to a comprehensive review [ 1149) . 

2.4.2. Helfrich's curvature Hamiltonian and its numerical implementation 

Another popular description consists of describing a membrane as an infinitely thin sheet that 
is characterized by a small number of elastic properties: spontaneous curvature and bending 
moduli. These coarse-grained parameters encode the architecture of the constituent lipids and 
the way they pack in the bilayer. In its simplest form the Hamiltonian is written as an expansion 




(50) 
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in the invariants that can be constructed from the local principal curvatures, Ci and C2 [ 

issiiissaisssi: 

7i= [d^A (-f + ^(Ci + C2- Cof + RC1C2) (51) 



2 

where Co is denoted the spontaneous curvature, and the integration is extended over the entire 
surface. The coefficient 7 is the lateral membrane tension, Kb and R, are the bending rigidity 
and saddle-splay modulus. This expansion is appropriate if the curvature is small compared to 
the characteristic curvature set by the inverse bilayer thickness or additional non-linear elastic 
effects. For homogeneous membranes of fixed topology, the integral over the Gaussian curvature, 
C1C2, is an invariant (Gauss-Bonnet theorem), and the term proportional to k contributes only 
a constant. The saddle-splay modulus is, however, important if the membrane topology changes, 
e.g., pore formation, fusion or fission. Membranes are often characterized as being tensionless, 
i.e., by 7 = 0. In this case, one can constrain the average membrane area to a fixed value, Aq, 
and utilize a term that is proportional to the square of the deviation of the membrane area. A, 
from it, 

AnA = Yi^-Aof (52) 

to describe the small area compressibility of the membrane. Such a model is able to describe 
successfully the large-length scale behavior of interfaces, membranes, and vesicles and their 
fluctuations. It has been the starting point for analytical calculations, such as the renor- 
malization of elastic constants by fluctuations of the local membrane position [ I334j . Fur- 
thermore, the model can very efficiently be studied by computer simulations. Much effort 
has been focused on exploring, for instance, the shape of vesicles and their fluctuations [ 

There also have been attempts to successively incorporate more local details into the descrip- 
tion by first applying the curvature Hamiltonian to each monolayer of a bilayer membrane and 
then to augment it by terms that account for the tilting of hydrocarbon tails [ I346| I347j . These 
techniques have been applied to estimate the free energy of localized, highly curved membrane 
structures as they occur at the edge of a bilayer [ I348j or in morphological changes of the bilayer 
structure, e.g., in transitions from the micellar to the inverted hexagonal phase [ I34(ij or in the 
fusion of membranes [ 349, 35(3 13511 1352j . The description of topological changes, e.g., pore 
formation, is difficult, however. 

One important simplification of the general curvature Hamiltonian H51() is achieved when 
overhangs of the membrane can be neglected. In this case the sheet can be described by a single- 
valued function, z = h{x,y). In this Monge representation the mean and Gaussian curvature, 
H and K, are given by 

_ C1 + C2 _ (1 + hl)hyy + (1 + hl)h^^ - 2h^hyKy 

2 " 2(l + /i2+/j2)3/2 ^^'^) 

and 

^ hxxhyy - hly 

where hx = dh/dx denote the derivative with respect to one of the two lateral coordinates. 

If the amplitude of fluctuations are small, then one can approximate 2H ~ h^x + hyy = Ah 
and Eq. (|^T|) takes the simple form (with Co = 0): 



n ^ I dxdy^/l + hl + hl{-f + -{Ahy) (55) 
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2int = j T^h exp 



(60) 



^ jA + J dxdy [j\^h\' + 2 (AM' j (56) 

Often this approximation is utilized to extract the interface tension and bending rigidity from 
computer simulations of small membrane patches [ 13531 11371 lllll 0511 11741 [T71 11821 1555) . To this 
end, one determines the local position of the membrane, h{x,y), from the simulation. Utilizing 
the Fourier spectrum of interface fluctuations, h(k), 

/i(x) = — r /t(k) expfik • x] = / d^k /i(k) exp[ik • x] (57) 

hik) = J d^x h{x) expl-ik ■ x] (58) 
one diagonalizes the Hamiltonian H56() 

= (^ / ^'"^ + f = ^ ^ + i^') 

which shows that the different Fourier modes of the local interface position are independent. 
The statistical mechanics of the interface position can be described by the partition function 

n[h] - 

'kBT_ 

where the functional integral / T>h sums over all local positions of the interface. Since the 
Hamiltonian is the sum of independent, harmonic degrees of freedom, /i(k), the Fourier modes are 
uncorrelated and Gaussian distributed around zero. Their variance is given by the equipartition 
theorem: 

^ [^yk^ + Kk') {\h{k)\^) = kBT (61) 

On short length scales, k > (j/ny^'^, the bending stiffness dominates and |/i(A;)p ~ k~^. The 
large scale fluctuations are controlled by the surface tension, 7, and obey \h{k)\'^ ~ k^"^. This 
method has successfully been employed to extract the tension and bending rigidity from simu- 
lations of interfaces in polymer blends [ 1159| 1356| 1357j , polymer-solvent interfaces [ 1317j liquid 
crystals [ ESHI ESU, and hpid bilayers [ ESSl HSU EH ESI CZll HIl CHI ESS]. To do so, one 
numerically determines the position, h, of the interface, or bilayer, in sub-columns centered 
around {xi,yj) and with lateral size A (block analysis [ 136U1 The distance between the 
grid points. A, is a compromise: if the discretization is too fine, the instantaneous composition 
profile along a column will strongly fluctuate. Composition fluctuations may affect the estimate 
of the interface position. If A is too large, however, the Fourier modes with large wavevector are 
averaged out and the spectrum is underestimated. This spectral damping is a function of the 
product of the spatial resolution A and the wavevector q. It only becomes irrelevant in the limit 
Aq <C 1. Cooke and Deserno [ 1175j explicitly calculated the strength of the spectral damping 
assuming that the smooth interface location is averaged inside each column. In this case, the 
spectrum of the grid-averaged fluctuations /i(k)gi.id and the original spectrum, /i(k) are related 
via 



|/i(k)g,id|' = \h{k)\' 



sin(AA;^/2) sin(AA;j^/2) 



{Ak^/2) {Akj2) J ^^^^ 

This discretization artifact becomes important when the fluctuation spectrum of very small 
membrane patches is studied by atomistic simulations or when the crossover between the tension- 
dominated and bending rigidity-dominated part of the spectrum occurs at rather large wave 
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Figure 5. Dynamic triangulation of fluid membranes. The bonds between vertices i and j has 
length lij, and the length of the dual tether, cjjj is indicated. The bond that connects vertices 
i and j is flipped. In order to conserve the triangular nature of the network each vertex must 
initially be connected to at least four neighbors. 



vectors (i.e., 7 > 0). A possible way to deal with these discretization artifacts is to divide out 
the a priori, /c-dependent spectral damping factor. 

The full Hamiltonian of Eq. (|51j) has to be used to describe large amplitude fluctuations 
of microemulsions and of closed vesicles. To this end one triangulates the surface. A typical 
snapshot is presented on the right hand side of Fig. ^ The vertices are connected by tethers 
that define the internal topology of the sheet. If the topology is fixed, the membrane will exhibit 
only elastic response that is characteristic, inter alia, of polymerized membranes. If one allows 
for changes of the internal topology in the course of the simulation (dynamic triangulation) [ 
13611 13361 1362j . one can mimic the in-plane fluidity. The vertices diffuse in this fluid membrane. 
A common strategy consists in cutting and re-attaching tethers between the four beads of two 
neighboring triangles. Such an elementary move is sketched in Fig. |S1 

The self-avoidance of the membrane is modeled by an excluded volume interaction between 
vertices at positions, Rj. The strength and range of the excluded volume interaction and the 
interactions along tethers can be chosen so as to avoid crossings. Let rij denote the surface 
normal of a triangle, i. The discretization of the mean curvature, H = n - AR at vertex i takes 
the form 

H, = -ni- J2 7^(R-*-Rj) (63) 

jdn.n.yi) 

where the sum runs over all tethers that are connected to the vertex i. The length of such 
a tether is lij, Oij = lij{cot @i + cot 02)/2 is the length of the corresponding tether in the dual 
lattice which is created from the intersections of the perpendicular bisectors of the bonds, ©i 
and G2 are the angles opposite to the link ij in the two triangles that border the link. The 
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quantity ai = \ J2jen.n.{i) ^ij^ij represents the area of the virtual dual cell. Using this expression 
and the fact that the local normal n is collinear to AR in three dimensional space, one obtains 



n = ^jd^AH^ = ^jd'A (AR) 



K X - 1 



51 -p(R-i-R-i 

J(in.n.{i) ^-^ 



(64) 
(65) 



for 7 = 0. We have omitted the integral over the Gaussian curvature which does not depend on 
the membrane configuration if the the topology is preserved. 

If one assumes that all triangles are equilateral, aij = lij/^f^, the equation above reduces to 



^ (riQ - rifif = \/3k ^ (1 - iIq • rvp) 



(66) 



where denotes the normal vector of the triangle a. In the general case of randomly trian- 
gulated surfaces, however, this latter expression suffers from deficiencies. Most intriguingly, use 
of the expression (|66() for a randomly triangulated sphere recovers the expected result, while a 
random triangulation of the cylinder in conjunction with Eq. (|66|) fails to yield the result of the 
continuum description [ 1363) . 

Generalizations of this class of models have be utilized to study freezing of vesicles [ I364| I365j 
or mixed membranes [ 13661 I367j . These models can also be employed in conjunction with a 
hydrodynamic description of the surrounding solvent [ I34U1 13681 13691 137Uj . 



3. An example of an integrated approach: fusion of membranes 

3.1. Motivation and open questions 

One example of a collective phenomena in membranes is the fusion of two apposing lipid 
bilayers. It is a basic ingredient in a multitude of important biological processes ranging from 
synaptic release, viral infection, endo- and exocytosis, trafficking within the cell, and fertilization 
[HlUniliniiSIlHlliniEni- The fusion process can be roughly divided into two steps [EOl: first 
the two membranes to be fused are brought into close apposition. Fusion peptides embedded in 
the membranes play an important role during this initial step. They ensure that only specific 
membranes are brought into close proximity with one another. One way to accomplish this is for a 
fusion protein, embedded in one of the membranes, to be is inserted into the apposing membrane, 
followed by a conformational change of the protein. This active mechanism imparts energy 
into the system. The specific role of proteins in the fusion process, their spatial arrangement 
and conformational changes, have attracted much interest for they are of great importance in 
regulating fusion [ 27, 371 . The second step consists of the fusion event itself in which the 
topology changes from two apposing bilayer of two vesicles to a fusion pore in a now-single 
vesicle. There is evidence that this second step is dictated by the amphiphilic nature of the 
bilayer constituents [ 13721 1261 1373j , for fusion occurs in very different systems ranging from 
tiny, highly curved, synaptic vesicles to whole cells. Moreover, sophisticated fusion peptides 
are not necessary to initiate fusion between laboratory vesicles. The simple depletion force 
that arises on the addition of water-soluble polymer (polyethylene glycol, PEG) to a vesicle 
solution [ 13741 13751 1376j . shear [ 1377j . or sonication [ 1377j serve equally well in inducing it. 
Another important piece of evidence comes from synthetic, polymer-based, membranes made 
of amphiphilic polymers [ 13781 11621 13791 138U1 13811 1382j . The behavior of these polymersomes 
resembles those of the much smaller and more fragile vesicles comprised of biological lipids and. 
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in the absence of proteins, includes processes like fusion and fission [ 13791 1577] . This all suggests 
that fusion is a universal collective phenomenon. Therefore coarse-grained models are well-suited 
to investigate its underlying mechanism. 

While many specific details are known about the first step of the fusion process, much less is 
known about the second. Generally, fusion is considered a "messy" process because of the drastic 
disruption of the bilayers' integrity. The time and length scales exclude a direct experimental 
observation of an individual fusion event. However even in the absence of direct information 
about the fusion intermediate, much can be inferred from a systematic variation of parameters 
(e.g., composition of mixed bilayers, or tension), and careful experimental techniques (e.g., 
electrophysiological measurements of membrane conductance). Some of the main experimental 
observations are: 

1. Lipids that favor a large negative spontaneous curvature of a monolayer, such as DOPE, 
increase the fusion rate [ [23 1383j . 

2. Increasing the tension of the apposed membranes results in an increased fusion rate [ I384j . 

3. During fusion, the lipids in the two apposing cis leaves mix [ 13751 1?7^ I385j . 

4. In some experiments, fusion is leaky. This implies that, correlated in space and time with 
the fusion event, there is a passage from the interior of the vesicle to the outside solution 
[ I386ll387ll388ll389ll39n] . 

5. Some experiments report on the transient mixing of lipids between the cis, or most closely 
apposed, leaf and the trans, or farther leaf, of the other membrane [ I37()|IH7K] . This process 
is much faster than the flip-flop tunneling of lipids from one leaf to the other in an intact 
membrane. 

The first three observations can be rationalized by the classical model of membrane fusion 
proposed by Kozlov and Markin in 1983 [ 139H 1392] . They conjectured a fusion pathway and 
calculated the free energy barrier utilizing an effective interface Hamiltonian. The monolayers 
of the membranes were modeled as thin elastic sheets, and the description was augmented by a 
free energy penalty for the packing difficulties that arise in the intermediate structures. 

In their model the proper fusion process starts with the formation of a stalk (see Fig. EJ, a 
rotationally symmetric connection between the two apposing cis monolayers. Once the stalk has 
formed, the two inner cis leaves retract from it leaving a transmembrane contact that consists 
of a small circular membrane patch built from the two outer trans leaves. This transmembrane 
contact can expand radially to form an extended hemifusion diaphragm. The rupture of this 
diaphragm creates a fusion pore. The expansion of the fusion pore completes the process. The 
classical fusion model is able to rationalize the first three observations. (1) Lipids that favor a 
large negative curvature of a monolayer tend to form an inverted hexagonal phase, and this non- 
lamellar morphology shares common local geometrical features with the stalk intermediate. (2) 
An increase of the lateral membrane tension, or free energy per unit area makes more favorable 
the decrease in membrane area which fusion brings about. (3) The stalk and the outer rim of 
the hemifusion diaphragm establish a connection between the two inner cis leaves along which 
lipids between the two cis monolayers can mix by diffusion. The last two experimental findings, 
however, cannot be explained. At no time during the classical fusion scenario is there a path 
between the inside of either vesicle and the outside solution. Also there is never any direct 
connection between the inner and outer monolayers of a bilayer membrane. 

The geometry of the stalk intermediate and the free energy penalty associated with the packing 
frustration of the hydrocarbon tails is an input into the classical model. Earlier calculations 
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Figure 6. Sketch of the classical fusion path. 

used a toroidal shape of the stalk, and assumed that the thickness of the monolayers remained 
constant as a void formed. An ad-hoc upper limit of the free energy costs of chain packing 
was obtained utilizing the macroscopic surface tension of hydrocarbons. The total free energy 
barrier associated with the formation of a stalk was estimated to be 200 ksT, a value much 
too large for fusion to occur in soft matter systems [ I393j . Subsequent refinements of the 
model utilizing a catenoid shape and including tilt and splay of the hydrocarbon tails solved 
this "energy crisis" [ 13491 1352j and yielded a significantly lower barrier of 30 to 40 ksT. The 
large variation of the estimated barrier with the assumptions of the model suggests that the 
effective Hamiltonian representing the monolayers as thin elastic sheets might not be able to 
accurately describe the highly bent structures of the intermediate. Models that retain the notion 
of amphiphilic molecules, that incorporate the packing of the molecular conformations in non- 
lamellar geometries, and that do not assume the structure of the fusion intermediates, are better 
suited to provide direct insights. 

3.2. Model and techniques 

By virtue of the experimentally observed universality of the proper fusion event, and due to 
the concomitant time and length scales that are unattainable with atomistically detailed models, 
coarse-grained models are well-suited to investigate the basic mechanism of membrane fusion. 
Ideally, such a coarse-grained model would combine the following characteristics: 

1. It describes the architecture of the amphiphilic molecules. The parameters of the model are 
directly related to experimentally measurable characteristics. The change of the molecular 
conformations and the associated loss of entropy in a non-planar environment can be 
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calculated. 

2. It can be used to observe whether the fusion event does, or does not, occur without prior 
assumptions about the pathway. 

3. It can be solved with a computational technique which permits one to simulate a number 
of independent fusion events. Note that in experiment, several different outcomes of an 
encounter between two vesicles can be observed, e.g., adhesion, complete fusion, or rupture. 

4. Its various parameters, such as the lipid architecture and thermodynamic conditions, can 
be explored, and the free energy of the fusion intermediates can be calculated. 

The first requirement is best met by a detailed model that mimics as much of the local structure 
of the bilayer and the surrounding solvent as possible. Such models, however, are computation- 
ally very expensive and do not permit the systematic exploration of the fusion process or the 
calculation of free energy barriers. 

In our own investigation, we have chosen the bond fluctuation model (cf. Sec. 12. 2[) . Coarse- 
grained models strike a balance between specificity of description and efficiency of computation. 
The bond fluctuation model, in which the amphiphiles are modeled as diblock copolymers on 
a cubic lattice and the solvent consists of homopolymers of identical length, is certainly one of 
the coarser, and therefore more computationally efficient, models. Although the model captures 
only the universal amphiphilic characteristics of the membrane components, it provides a rea- 
sonable description of bilayer properties (cf. Tab.P). This efficient model allowed us to study 
rather large systems that contain several thousand amphiphiles, and to observe thirty-two fusion 
events for one set of parameters. Besides computational efficiency, the model has two additional 
advantages: (1) Much is known about the properties of the model, e.g., interface tension, bi- 
layer compressibility, phase behavior, spectrum of interface fluctuations, etc. (2) The model can 
be quantitatively compared to the standard Gaussian chain model of the SCF theory without 
any adjustable parameter. We use Monte Carlo simulations [ [5T1 1161] to provide an unbiased 
insight into the fusion pathway. Once this is attained, we perform extensive SCF calculations 
[ 131 H I312j of the same model in order to obtain the free energies of intermediates over a wide 
range of parameters. 

3.3. MC simulation 

3.3.1. Separation of time scales 

Our Monte Carlo simulations are performed in the canonical ensemble [ 1511 1161] . The molec- 
ular conformations are updated by local segment displacements and slithering-snake-like move- 
ments. These movements conserve the local densities and thus lead to a diffusive behavior on 
large length scales. Moreover, the molecules cannot cross each other during their diffusive mo- 
tion. In this sense we have a slightly more realistic time evolution on local length scales than 
in dissipative particle dynamics simulations [ I142j . but Monte Carlo simulations cannot include 
hydrodynamic flow, which might become important on large length scales. We count one at- 
tempted local displacement per segment and three slithering-snake-like attempts per molecule as 
four Monte Carlo steps (MCS). This scheme relaxes the molecular conformation rather efficiently 
[ 1394] . The time scale of the MC simulations can be compared to experiments by matching the 
self-diffusion coefficient of the lipids in a single bilayer (see below). At any rate, we do not 
expect the time sequence to differ qualitatively from that of a simulation with a more realistic 
dynamics on time scales much larger than a single Monte Carlo step. Most importantly, fusion 
is thought to be an activated process. Therefore the rate of fusion is dominated by free energy 
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barriers encountered along the fusion pathway, while the details of the dynamics only set the 
absolute time scale. 




Figure 7. Total thickness of the bilayer membrane, d, measured in units of the radius of gyration, 
Rg ~ 7u, vs. the exchange chemical potential Afi between amphiphiles and solvent molecules. 
Membrane tension, 7, as a function of exchange potential, A/U, is shown on the right scale. 



We begin our simulation by preassembling flat, tense, bilayers. The tension is dictated by 
the areal density of amphiphiles, i.e., the thickness of the membrane (cf. Fig. EJ. In the lattice 
MC simulations we cannot simulate at constant pressure or surface tension (which is routinely 
used in off-lattice models [ KS95( KS96j ) because the size of the simulation box cannot be changed 
continuously. The lattice also prevents us from measuring the pressure or tension via the virial 
expression (cf. Eq. (|33j) ) or the anisotropy of the virial of the forces. The thickness of a tensionless 
membrane, however, can be measured by simulating a bilayer patch that spans the periodic 
boundary conditions of the simulation cell in only one direction but exhibits a free edge in the 
other direction. The area of the bilayer then will adjust so as to establish zero tension and 
the thickness can readily be measured in the center of the bilayer patch. Such a tensionless 
configuration is shown in the middle of Fig. ^ 

The dependence of the tension on the exchange potential, A/i, between amphiphiles and sol- 
vent can be obtained by simulations in the grand-canonical ensemble where the density is con- 
stant but, in addition to the MC moves that renew the molecules' conformations, one "mutates" 
amphiphiles into solvent molecules and vice versa [ I,'ffl7| 0398 . These moves change the composi- 
tion of the system while the thermodynamically conjugate variable, the exchange potential A^, 
is controlled. 

The excess free energy per unit area, in the thermodynamic limit of infinite area, defines the 
lateral membrane tension 

A^oo A. 

5nm(T,Afi,A) = n^{T,Afi,V,A)-no{T,Af,,V) (67) 

where A denotes the area of the membrane, and 0,m and the grand canonical free energy of 
the system with and without membrane, respectively. In an incompressible system the tension 7 
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can be related to the temperature and chemical potential by means of the Gibbs-Duhem equation 

d7(r, Afi) = -6sdT- 6aad{Afi) , (68) 

where 5s is the excess entropy per unit area, and 6(Ta is the excess number of amphiphilic 
molecules per unit area. This relation can be exploited to integrate the thickness of the bilayer 
d = aaN/ p with respect to the exchange potential Afi at constant temperature and calculate the 
change of the membrane tension [ I353j . Using the tensionless state as staring point we obtain 
the relation between the membrane thickness, d and the tension, 7 as shown in Fig. 
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Figure 8. Evolution of internal energy in fusion simulations. The curves correspond to dif- 
ferent initial bilayer separations A = 4n (squares) and A = lOn (circles) and bilayer tensions, 
7/7int = 0.75 and 0.375, respectively. To reduce fluctuations, the data are averaged over all 32 
configurations at equal time and additionally over small time windows. The large negative value 
of the energy mirrors the attractive interactions in the solvent. The inset shows the sample- 
to-sample energy fluctuations as a function of time (for the system with high tension and large 
membrane separation). Large fluctuations identify the onset of fusion. 



After the bilayer has assembled without defects in a thin film geometry, we stack two bilayers 
on top of each other with a distance A between them. The system of two apposing bilayers is 
embedded into a three dimensional simulation cell with periodic boundary conditions in all three 
directions. Configuration bias MC moves [ 13991 11001 l4Ulj are utilized to fill the remaining volume 
with homopolymer solvent. This starting configuration of flat bilayers mimics the approach of 
two vesicles whose radii of curvature are much larger than the patch of membrane needed for 
fusion. Then we let the system evolve and observe the fusion process in the simulation. 

Because fusion is a kinetic, irreversible process, the starting condition might have a pronounced 
influence on the outcome. In Fig. IS] we present the internal energy of the system as a function 
of MC steps on a logarithmic scale. Immediately after the assembly of the two bilayer system, 
one observes a rapid decrease of the energy. This corresponds to the relaxation of the bilayer 
structure after the insertion of the solvent. After r/j ~ 25000 MCS (~ 9ns, see below), however. 
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the local internal structure of the bilayer has equilibrated and the energy starts to rise very 
slowly. This increase is compatible with a logarithmic growth law and stems from thermally 
excited interface fluctuations that increase the area of the hydrophilic-hydrophobic interface. At 
the end of this stage, stalks and holes are formed in the simulations. (Here and in the following, 
we denote a pore across a single bilayer as a "hole" in contrast to the "fusion pore" that spans both 
bilayers.) Finally, around Tp ^ 8 - 10® MCS, the energy suddenly drops indicating the formation 
of a fusion pore that reduces the total bilayer area and, therefore, the stored tension. The 
important point is that there is a difference of two orders of magnitudes between the relaxation 
time of the local structure of an individual bilayer, and the time scale on which fusion occurs. 
Due to this separation of time scales between initial relaxation and fusion, we do not expect 
the preparation of the starting configuration to affect the irreversible fusion process. Similarly 
we do not expect our results to depend on our particular choice of relaxation moves, as other 
choices would also lead to relaxation of the bilayers which takes place on a much shorter time 
scale than does fusion. 

The inset of Fig. |H1 shows the fluctuations in the energy, i.e., the fluctuations between the 
thirty-two different runs at equal time. Strong fluctuations indicate energy differences between 
the independent runs. The peak at around Tp ~ 8T0® MCS (for the system with high tension and 
large separation) indicates that some systems have already formed a fusion pore, and therefore 
have a lower energy, while other systems have only stalks and holes, and therefore have a higher 
energy. The width of the peak provides an estimate for the spread of the time during which 
a fusion pore appears. Specifically, for the system with the lower tension and large initial 
separation, at time 1.26 • lO'^ MCS we have observed the formation of stalks in 23 systems while 
9 runs have already formed a fusion pore. At time 2.52 • 10'' MCS, 20 runs have successfully 
completed fusion, in 3 configurations only stalks have formed and in 9 runs stalks and holes 
have appeared but the fusion has not been successfully completed. The figure shows data for 
two different bilayer tensions and two different distances between the bilayers. Increasing the 
bilayer tension, and reducing the bilayer distances accelerates the fusion process. Moreover, the 
slower the fusion process, the clearer is the transition state we observe. Fusion occurs around 
8 • 10® MCS and 2 • lO'^ MCS for the systems at large separation and with high and low tension 
respectively. 

The lateral diffusion constant of a lipid in a single bilayer is I? ~ 10~^?/^/MCS. If we identify 
the thickness dc = 21u of hydrophobic core in the tensionless state with 3nm and utilize a 
typical value for the self-diffusion coefficient of lipids in bilayer membranes, D « 6 • 10®nm^/s, 
we estimate that one MCS corresponds to 0.36 ps. Thus, the time scale on which we observe 
fusion is 3//s and 7/is for the tense and less tense bilayers, respectively. This estimate is about 
an order of magnitude larger than what is observed in the DPD simulations of Shillcock and 
Lipowsky [ ,402^ and about 2 orders of magnitude slower than the fusion process in the Lennard- 
Jones-type model of Smeijers et al [ I403| . In the latter study, however, the time scale was not 
estimated from the comparison of the lipid self-diffusion coefficient in a bilayer but directly from 
the time of the MD simulations without further rescaling which is likely to underestimate the 
duration of the fusion process [ 404 . 

Unless noted otherwise, we will discuss in the following the data for the larger tension and 
larger initial bilayer separation. The data for the lower tension are very similar. 

An important property of the bilayer which we prepare is its tension, or free energy per unit 
area. Thermodynamically, a tense bilayer is only metastable, and will eventually rupture. Hole 
formation reduces the bilayer's area and lead to a stable, tensionless state. Hole formation is 
an activated process and, indeed, we do observe the rupture of isolated tense bilayers in very 
long simulation runs. Since vesicles have to be stable for long times in order to fulfill their 
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carrier and enclosing function, it is reasonable to require that the time scale of hole formation 
and rupture, th, of an isolated membrane be much larger than the times scale of fusion of two 
apposed bilayers. Yet in order to undergo fusion, just such long-lived holes must occur at some 
point along the fusion pathway. It would seem that vesicles could either be stable, or they could 
undergo fusion, but not both. The question that immediately arises is how membranes actually 
manage to exhibit these two conflicting properties. We return to this question in Sec. |3 
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Figure 9. Area of holes vs. time in the sys- 
tem of two apposed bilayers (gray for one 
bilayer and black for the other on the bot- 
tom panel) and in an isolated bilayer (top 
panel). From Ref. [ I51j 



In Fig. 1^ we present time traces of the areal fraction of holes in an isolated tense bilayer and 
in the two apposed bilayer system. In the isolated bilayer, small holes form and close but the 
bilayer remains stable on the time scale of the fusion process. In the two apposed bilayer system, 
however, larger holes form more readily and their occurrence is correlated with the fusion event. 

This additional MC data further supports the observation that in our coarse-grained model, 
there is a clear separation of time scales, r/j <C t^t <^ th, between the local relaxation time of 
bilayers, tr, the time scale, Tp, to nucleate a fusion pore, and the lifetime of an isolated tense 
bilayer, th- 

3.3.2. Observed fusion pathways 

During the initial stage of simulations, interface fluctuations of the initially flat bilayers are 
thermally excited and the bilayers collide with one another. Sometimes these collisions give rise 
to small local interconnections. For the most part, these contacts are fleeting. Occasionally we 
observe sufficient rearrangement of the amphiphiles in each bilayer to form a configuration, the 
stalk, which connects the two bilayers. Once such a stalk has formed, it is rather stable on the 
time scale Tp. Density profiles of the hydrophilic and hydrophobic parts of the amphiphiles in 
the presence of the stalk, obtained by averaging over configurations, are shown in Fig. IIUI The 
dimples in the membranes at each end of the stalk axis are notable. What can barely be seen is 
a slight thinning of each bilayer a short distance from the axis of the stalk. 

Under the specific conditions of our simulations we believe that stalks are only metastable. 
First, the observed stalks are isolated and their occurrence goes along with a rise of the internal 
energy. Second, occasionally we observe that stalks disappear without proceeding further to a 
fusion pore. This behavior indicates that the stalk represents a local (metastable) minimum 
along the fusion pathway. Once a stalk has formed, we observe that it elongates asymmetrically, 
and moves around in a worm-like manner. Evidently its elongation does not cost a great deal 
of free energy. 
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Figure 10. Density distribution of segments 
in the stalk, averaged over all 32 simula- 
tion runs. At each point only the major- 
ity component is shown: solvent as white, 
hydrophobic and hydrophilic segments of 
amphiphiles as black and gray respectively. 
From Ref. [EI] 



After stalks are formed, the rate of formation of holes in either of the two bilayers increases. 
Stalk and hole formation are not only correlated in time but also in space: Holes form close to 
the stalks. This can be quantified by the hole-stalk correlation function: 

J2r,,rJ{\rs-rh\-r)Pshirs,rh) 

air) = 77] , ^ (69) 

Er^r;, 0{\rs-rh\ -r) 

where Ps/j(rs,r/j) is the joint probability that the lateral position is part of a stalk and 
is part of a hole, and 5{r) is the Dirac delta function. The value of g{r) at large distances r is 
proportional to the product of the areal fraction of holes and stalks. This correlation function is 
shown in Fig. 1121 The scale of g{r) increases with time indicating the simultaneous formation of 
stalks and holes. The figure shows that the correlation peaks at a distance of about 2/3 bilayer 
thickness, and rapidly decays at larger distances. 

Now that a stalk connecting the two bilayers has appeared as well as a hole in one of the bi- 
layers, fusion proceeds along one of two closely related paths. These are depicted in Fig. ^2 The 
snapshots are taken from two representative simulation run. Configurations are numbered by the 
time (in multiples of 25 000 MCS) at which it was observed. Three-dimensional plots of the local 
density are shown - the hydrophobic core is depicted as dark gray, the hydrophilic-hydrophobic 
interface (defined as a surface on which densities of hydrophilic and hydrophobic segments are 
equal) is light gray. Hydrophilic segments are not shown for clarity. Each configuration is shown 
from four different viewpoints. Top- and bottom- left sub-panels have been generated by cut- 
ting the system along the middle plane parallel to the plane of the bilayer. The top and bottom 
halves are viewed in the positive (up) and negative (down) direction correspondingly. In these 
panels one sees the cross-section of connections between bilayers (hydrophobic core of stalks) 
as dark regions. Holes in the individual bilayers appear as white regions in one of the panels. 
Top- and bottom- right sub-panels are side views with cuts made by planes perpendicular to the 
bilayer. The grid spacing is 20n dc (see Tab.^). 

1. A hole appears in one bilayer and the stalk completely surrounds it rather rapidly. The 
resulting configuration looks very much like a hemifusion diaphragm which has been sug- 
gested by many authors as an intermediate stage in fusion [ .392. ,,405.. .393^ . However, this 
diaphragm is quite different from the hemifusion diaphragm of the classical scenario that 
consists of two trans monolayers of the fusing membranes. In contrast, the diaphragm we 




Figure 11. Two observed pathways of fusion process. The snapshots are taken from two repre- 
sentative simulation runs. Each configuration is numbered by the time (in multiples of 25, 000 
MCS) at which it was observed. For discussion of the mechanism see text. From Ref. [ I13j 
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Figure 12. The hole-stalk correlation function, g{r), at early times. From Ref. [I51j. 



observe in our simulations is made of one of the pre-existing bilayers; i.e., it is comprised 
of cis and trans leaves of the bilayer that did not form a hole. The appearance of a pore 
in this diaphragm and its expansion completes the formation of the fusion pore. 

2. A hole appears in one bilayer and, before the stalk completely surrounds it, a second 
hole appears in the other bilayer. The stalk continues to surround them both, and align 
them both during this process. Eventually, the stalk completely encircles both holes and 
a complete and tight fusion pore is formed. This path is slower than the previous one and 
the holes tend to expand more during this process. 

Once the fusion pore has formed, it expands because it reduces the bilayer area and thereby 
relieves the tension. Note that in the canonical ensemble, the total tensionless area is fixed from 
the beginning and the growth of the fusion pore ends when the pore reaches its optimum size 
determined by the finite size of our simulation cell. Very similar finite-size related effects have 
been studied for the hole formation in canonical simulations of isolated bilayers [ l4U6j . 

Occasionally, more than one stalk forms in the simulated bilayer patch. An example of con- 
figurations with multiple stalks is shown in Fig. 1131 The interactions between stalks have been 
considered by Lukatsky and Frenkel within an effective interface model [ 14071 . They argue that 
membrane- mediated interactions lead to an attraction of stalks, and that the collective cluster- 
ing of stalks, in turn, aids the fusion process. This is compatible with the snapshots presented 
in Fig. 1131 where there is apparently a higher probably of finding two stalks close to each other. 

The observed pathway of fusion differs markedly from the classical hypothesis in the following 
important aspects: First, the fusion intermediates we observe break the axial symmetry which 
has been assumed in previous calculations. Second, holes in the individual bilayers play an 
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Figure 13. Snapshot from a simulation 
with small membrane tension and large bi- 
layer distance (cf. Fig. l21() showing multiple 
stalks and a hole. The configurations are 
depicted as in Fig. 
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Figure 14. (left) Time dependence of the area of holes (in single bilayers) and the fusion pore in 
one (out of 32) representative simulation run. Note the different scales for hole and pore areas. 
From Ref. [ |^ (right) Electrophysiological experiment on influenza hemagglutinin-mediated 
fusion of HAb2 and red blood cells, upper panel: equivalent electrical circuit, lower panel: 
Fusion experiment showing leakage temporarily correlated with fusion. The pore conductance 
Gp = 1/Rp marks the opening (arrow), flickering and growth of the fusion pore. Gdc is the 
conductance of the HAb2 cell when the fusion pore is closed and is the sum of the conductances 
of the HAb2 and the red blood cell when the fusion pore is open. It is a measure for the leakiness 
of the fusion event. Ten out of twelve experiments showed leakage. From Ref. [ I390j 
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important role in the fusion mechanism. On the one hand, holes give rise to some degree of 
transient mixing of the amphiphiles in the trans and cis leaves of the same membrane that is 
correlated with the fusion event. Lipids can switch between the two leaves by diffusion along 
the rim of these holes. On the other hand, the formation of holes in the individual membranes 
implies a transient content leakage that is correlated in space and time with the fusion event. 
Such leakage has been observed in recent experiments [ I388| 13^9(1] . In Fig. El we present the areal 
fraction of holes in the individual membranes and the size of the fusion pore from one simulation 
run. One clearly observes that hole formation precedes the fusion pore. The extent of leakiness 
that can be observed in experiments depends on the substance that passes through the holes 
of the vesicles and the fraction of the rim of a pore that is sealed by the stalk. If the stalk is 
very elongated and covers a substantial fraction of the incipient hole then leakage will be very 
small. In the same figure, we present electrophysiological experiments by Frolov and co-workers 
[OHHlinnni that show the conductance between two fusing vesicles and the conductance between 
the individual vesicles and the solution. The former quantifies the size of the fusion pore, while 
the latter indicates the area of holes in the individual membranes. Consistent with the Monte 
Carlo simulations there is no connectivity between the inside of the vesicle and the outside 
solution before the fusion event, i.e., the vesicles are tight. Just before a current between the 
vesicles is observed, however, the experiments reveal a substantial leakage current. Although 
the simulation and experiments deal with quite different systems, they both observe leakage in 
contrast to the classical hemifusion hypothesis. This exemplifies the insights into the mechanisms 
of collective phenomena that one can gain from simulations of coarse-grained membrane models. 

3.3.3. Comparison to other coarse-grained models 

While this qualitative agreement with experiment is very encouraging, the coarse-grained 
simulation model differs in many aspects from experiment. In Sec. 12.21 we demonstrated that 
our model can reproduce many properties of bilayer membranes, and in the previous section we 
have argued the the time scale of fusion is clearly separated from the relaxation time of local 
bilayer properties. Therefore we do not expect that including details of the lipid architecture or 
details of the diffusive dynamics of the MC simulations will qualitatively alter our conclusions. 

Nevertheless to gauge the universality of the observation of this specific coarse-grained model, 
it is very important to relate the findings to results of alternative coarse-grained models. Since 
the various models include different details, they emphasize different aspects of the fusion process. 
The simulation studies of fusion differ both in the geometry of the initial state - two apposed 
planar bilayers, or a planar bilayer in contact with a vesicle or two small vesicles ~ as well as in 
the representation of the amphiphiles and the solvent. 

Chanturiya et al. [ I408j investigated the role of tension on fusion within a two-dimensional 
model of a lipid bilayer. This limitation excluded the possibility of observing complex fusion 
intermediates. 

Noguchi and Takasu [ ll4Uj studied the fusion of two small bilayer vesicles using a solvent- 
free model [ I171j (see Sec. 12. 2|) . The amphiphiles were modeled as a rigid rod comprising one 
hydrophilic and two hydrophobic beads. Density dependent potentials were tuned to bring 
about the self-assembly in the absence of solvent. One important difference to the fusion of 
two planar bilayer membranes consists in the rather small contact zone of the two vesicles and 
the high curvature of the bilayer membranes. They observed two distinct fusion mechanisms: 
One mechanism resembled the classical hemifusion mechanism starting from the formation of a 
stalk and subsequent transmembrane contact. The transmembrane contact, however, did not 
significantly expand but remained confined to the rather small contact area of the two vesicles. 
A pore in the transmembrane contact completed the fusion. The other mechanism they saw was 
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Figure 15. Sequential snapshots of the 
fusion of two vesicles within the solvent- 
free model of Noguchi and Takasu [ I171j . 
Hydrophilic and hydrophobic segments are 
represented by spheres and cylinders, re- 
spectively. Panel (a) presents a snapshot of 
the initial stalk between the two apposed 
vesicles. Panels (b), (c), (d) show cuts per- 
pendicular to the plane of contact while 
the images (b'), (e), (f), and (c') depict 
cuts parallel to the plane of contact where 
the elongation of the stalk is clearly visible. 
From Ref. [HHJ. 
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the stalk-hole mechanism which we had observed. Their work was completely independent of 
ours, and was carried out at about the same time. In particular, they observed the elongation 
of the stalk along the edge of the contact zone between the two vesicles. These observations are 
visualized in Fig. [T5l 

The self-assembly of vesicles from a disordered solution and their fusion has been studied using 
MESODYN simulations [ l4fl9|I^Tf1] of diblock copolymers [ I411j . The model is similar to the one 
used in Sec. 13.41 but the use of a very coarse discretization of the molecular architecture permits 
the study of large system sizes in three dimensions. Interestingly, these calculations also observe 
the formation of stalks between two apposed vesicles and the formation of holes in the vicinity 
of stalks in both bilayers (stalk-hole mechanism 2). The subsequent encircling of both holes by 
the stalk completes the fusion process. 



ft) 




Figure 16. Sketch of stalk formation due 
to association of splayed lipids suggested by 
coarse-grained simulations, (a) Two vesi- 
cles are brought together and (b) a flat con- 
tact forms where, at the edge, the area per 
lipid in the outer leaflet increases as the 
membrane is strained, (c) Lipids tilt at 
the contact and (d) the hydrophobic tails 
of some amphiphiles begin to splay, (e) 
Splayed molecules then associate by their 
tails to form a new hydrophobic core, (f) 
which expands as the tails extend to form 
a classical stalk-like structure. From Ref. [ 

En]. 



Stevens, Hoh and Woolf [ I412j also studied the fusion of two small vesicles that had been 
pushed together via an external force. The double-tailed amphiphiles were described by a bead- 
spring model of eleven particles, and the solvent as consisting of single particles. Their simulation 
showed a single pathway to fusion, the stalk-hole mechanism. As in Noguchi and Takasu's second 
mechanism, they found a highly asymmetric expansion of the stalk along the edge of the contact 
zone. Most notably, the simulations provided a detailed description of stalk formation. Since 
the hydrocarbon tails were modeled as semi- flexible chains, they were found to tilt at the rim of 
the contact zone between the vesicles. The model also afforded the possibility of double-tailed 
lipids bridging between the cis-layers of the two apposed vesicles. The authors argue that this 
effect is important for nucleating the initial stalk. The details of the stalk formation inferred 
from the simulations are sketched in Fig. 1161 

Utilizing a similar coarse-grained bead-spring model [ I413j . Smeijers et al [ I4f)3j also provide 
a detailed description of stalk formation and emphasize the role of fluctuations that give rise to 
microscopic hydrophobic patches of the bilayers. In contrast to the simulations of Stevens et al 
[ I412j . however, stalks do not necessarily form at the edge of the contact zone. One important 
difference between Smeijers' simulations and others is the absence of tension in the vesicles that 
have been prepared by spontaneous self-assembly from a disordered solution. For the small 
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Figure 17. Fusion of two tensionless vesicles observed in a coarse-grained bead-spring model 
[ I413j . Shown is a cross-section along the vesicle- vesicle axis. Solvent is shown as dark blue 
spheres, and the solvent particles that enter the bottom vesicle highlighted light blue, (a) Some 
solvent enters the vesicle and initializes a hole in the lower vesicle (26.5 ns). (b) The hole is 
evident (26.8 ns). (c and d) As the stalk encircles the hole and the last solvent particles enter 
the vesicle, a hemifusion diaphragm is formed by the two monolayers of the upper vesicle (27.0 
ns and 27.2 ns). From Ref. [ i403) 

vesicles considered in the simulations, the curvature suffices to induce fusion. For two planar 
bilayers they observed the formation of stalks but no fusion [ !4n4j . In these simulations the stalk- 
hole mechanism was observed (see Fig. I17p significantly more often than the radial expansion of 
the stalk. 

Even more details of the lipid architecture are incorporated in the systematically coarse- 
grained representation of Marrink and Mark [ who studied fusion of two very small vesicles 
[ I414j . This model corroborates the tilting of the amphiphilic tails in the stalk. Marrink and 
Mark observe two pathways: the classical hemifusion mechanism (pathway I) and the stalk-hole 
mechanism, with an elongated stalk (pathway II). Importantly, transient pores in one of the 
bilayers were observed in the vicinity of the stalk in the second mechanism. The two observed 
pathways are depicted in Fig. ITHl 

The simulation model of Shillcock and Lipowsky takes a step towards an even more coarse- 
grained approach utilizing a DPD-model [ I142j . This allows them to study a large number 
of fusion events between a tense planar bilayer and a tense vesicle [ I4U2| . The larger degree 
of coarse-graining also allows them to systematically vary the tension of the two membranes. 
While they do not discuss the detailed mechanism of the fusion, they find that successful fusion 
is limited to a narrow range of rather high tensions. If the tension is too high, however, the 
membranes rupture instead of fusing. If the tension is below some threshold, the two apposed 
membranes do not fuse but rather adhere on the time scale of the simulations. This strong 
tendency to adhere has not been reported in other simulation models. Thus it does not appear 
to be a universal characteristic, and it remains unclear which specific property of the model is 
responsible for this feature. 

Utilizing a DPD-model, Li and Liu investigated the structure of the fusion intermediates [ I415j . 
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Figure 18. Comparison of the transition from stalk 
to opening of the fusion pore in both pathways ob- 
served for the fusion of mixed DPPC/DPPE vesi- 
cles. Slabs perpendicular to the fusion axis are 
shown, cutting through the stalk or hemifusion di- 
aphragm. Lipid headgroups are represented by large 
spheres. Different colors distinguish between lipids 
in the inner (yellow/silver) and outer monolayer 
(brown/black) and between the two vesicles (brown 
and yellow vs. black and silver). Orange spheres de- 
note the ethanolamine site of PE, blue spheres denote 
exterior water, purple spheres interior water. Note 
the differences in stalk structure (bent in pathway II) 
and in the composition of the HD (mixed in pathway 
I, almost entirely from a single vesicle in pathway II). 
From Ref. [ I414j (suppl. information). 



They found an asymmetric elongation of stalks similar to our observation. For rather symmetric 
amphiphiles, however, the elongated stalks expanded into an axially symmetric transmembrane 
contact, while for more asymmetric lipids, holes formed in the individual membranes in four out 
of five simulation runs. 

Taken together, the different coarse-grained simulation studies provide a consistent picture of 
the microscopic details of membrane fusion. In particular, the stalk-hole mechanism appears to 
be a viable alternative to the classical hemifusion mechanism. Although the different lattice, 
bead-spring, and soft DPD models differ substantially in their microscopic interactions, and 
in the fusion geometry, planar bilayer or highly curved vesicle, they all observe the non-axially 
symmetric elongation of the stalk and transient holes in the individual membranes in the vicinity 
of a stalk. This is also compatible with the experimental observation of transient leakage in some 
experiments. 

It is difficult to determine from the above results the conditions under which the classical 
fusion mechanism or the stalk-hole mechanism is the favored one. Field-theoretic methods, 
however, are well suited to explore model parameters and to answer this, and other, questions. 

3.4. SCF calculations 

One advantage of our coarse-grained model is the fact that it can be mapped onto the standard 
Gaussian chain model that is routinely utilized in SCF calculations [ I288|l285] . The length scale. 
Re, and the incompatibility xA^, can be extracted from the simulations. Further, simulation 
results of the bond fluctuation model have quantitatively been compared to SCF calculations 
for many properties [ 416 . 4T71 14181 UTOl fl58 . 356, 25^. The degree of the quantitative agreement 
without any adjustable parameter can be gauged from Fig. [3 

The SCF theory approach allows us to calculate the free energy of axially symmetric inter- 
mediates utilizing the radius of the structures as a reaction coordinate. First, we examine the 
dependence of the intermediates along the classical hemifusion path as a function of the bilayer 
tension and lipid architecture. Then, we proceed to estimate the free energy of the alterna- 
tive fusion intermediates observed in the simulations by patching together radially symmetric 
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structures. These SCF calculations go beyond the phenomenological approaches that utilize ef- 
fective interface Hamiltonians in two crucial aspects: (1) They retain the notion of amphiphilic 
molecules and calculate the changes of the molecular conformations in the complex, spatially 
inhomogeneous environment. No assumptions about the chain packing or stretching have to 
be made. (2) Only the radius (reaction coordinate) and the topology of the intermediate is 
dictated. The detailed geometry of the intermediate structure is optimized as to minimize the 
free energy that comprises contributions both from the interface between hydrophobic and hy- 
drophilic components, and from the loss of configurational entropy due to changes of the chain 
conformations. 

3.4.1. Axially symmetric configurations along the classical fusion pathway 

The free energy of an axially symmetric stalk that connects the two apposed bilayer membranes 
as a function of its radius is shown in Fig. 1191 The free energy profiles exhibit two maxima - the 
first one, ^o, corresponds to the formation of an initial connection between the bilayers (stalk) 
and the second one, 5*2, corresponds to the expansion of the stalk to a hemifusion diaphragm. 
These two maxima are separated by a minimum, Si, that marks the radius of a metastable 
stalk. The radius of this metastable configuration is set by the bilayer thickness. Density 
plots of the majority component of the three extremal states are shown in Fig. 1191 While the 
phenomenological calculations focused much effort on calculating the free energy cost of forming 
a stalk, So, the SCF calculations show that the main barrier along the path towards the expanded 
hemifusion diaphragm is not the formation of the initial stalk, Sq, but it is determined by the 
expansion, 52, of the metastable stalk to the diaphragm. 

The free energies of the metastable stalk. Si, and the saddle point, 5*2 along the path towards 
the hemifusion diaphragm are shown in Fig. [20] as a function of the amphiphilic architecture, 
/, and the bilayer tension, 7. The free energy barrier, S2, associated with the stalk expansion 
strongly decreases with tension. Within the classical model, this explains why the fusion rate 
increases with tension. The free energy of the metastable stalk hardly depends on tension 
but it decreases substantially as we lower the fraction / of hydrophilic segments and thereby 
decrease the spontaneous curvature of a monolayer to more negative values. For very asymmetric 
amphiphiles, those that form the inverted hexagonal phase in the bulk, the free energy of a stalk 
actually becomes negative. In this case, many stalks are expected to form and to condense on an 
hexagonal lattice thereby creating a "stalk phase". Such hexagonally perforated lamellar phases 
have been observed both in diblock copolymers [ I42U1 and in lipid/water mixtures [ I422j . 
Making the amphiphiles even more asymmetric, we observe that stalks spontaneously elongate 
and the inverted hexagonal phase forms. On the other hand, if we make the amphiphiles more 
symmetric, the metastable stalk configuration eventually disappears. In this case, the fusion rate 
is not determined by the free energy difference between the saddle point S2 and the metastable 
stalk Si but by the large barrier S2 only. Thus, the absence of a metastable stalk will strongly 
suppress fusion. 

Both in the classical mechanism, as well as in the alternate stalk-hole fusion mechanism 
observed in the simulations, stalks play a pivotal role. The limits of metastability of stalks 
explored by the SCF calculations are compiled in Fig. 1211 Most notably, the free energy calcula- 
tions demonstrate that fusion is restricted to a rather narrow window of amphiphilic architecture 
characterized by the ratio of spontaneous monolayer curvature and bilayer thickness. This nar- 
row range of spontaneous curvatures extracted from a coarse-grained model is bound by the 
spontaneous curvatures of two relevant lipids, DOPE and DOPC. 

The rupture of the hemifusion diaphragm completes the fusion process. The free energy of 
a fusion pore in tensionless bilayers as a function of its radius, R, is shown in Fig. For 
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Figure 19. (a) The free energy, F, of the stalk-Hke structure connecting bilayers of fixed, zero, 
tension is shown for several different values of the amphiphile's hydrophilic fraction /. In the 
inset we identify the metastable stalk. Si, the transition state, 5*0, between the system with no 
stalk at all and with this metastable stalk, and the transition state, 52, between the metastable 
stalk and a hemifusion diaphragm. The architectural parameter is / = 0.30 for this inset. No 
stable stalk solutions were found for / = 0.45 in the region shown with dashed lines. They 
were unstable to pore formation. Profiles showing the majority component of the two barriers 
and 5*2, and the metastable stalk, 52, are shown below the main panel, (b) The free energy 
of the expanding stalk-like structure connecting bilayers of amphiphiles with fixed architectural 
parameter / = 0.35 is shown for several different bilayer tensions. These tensions, 'j/^mt, are 
shown next to each curve. From Ref. [ I311j . 
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Figure 20. The free energy, F, of 
the metastable stalks Si (dashed hnes) 
and the transition states 5*2 (full lines) 
along the path towards the hemifusion di- 
aphragm are shown clS Si function of the 
tension for different architectures / = 
0.25, 0.30, 0.35, 0.40. Notice that there 
is no metastable stalk (solution Si) for / = 
0.4 at zero tension. From Ref. [ I311j . 




Figure 21. A "phase diagram" of the fusion process in the hydrophilic fraction-tension, (/, 7), 
plane. Circles show points at which simulations were performed by us. Successful fusion only can 
occur within the unshaded region. As the tension, 7, decreases to zero, the barrier to expansion 
of the pore increases without limit as does the time for fusion. As the right-hand boundary is 
approached, the stalk loses its metastability causing fusion to be extremely slow. As the left- 
hand boundary is approached, the boundaries to fusion are reduced, as is the time for fusion, 
but the process is eventually pre-empted due to the stability either of radial stalks, forming the 
stalk phase, or linear stalks, forming the inverted hexagonal phase. From Ref. [ I311j . 
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Figure 22. The free energies, -F, of a fu- 
sion pore (solid lines) and of a stalk (dashed 
lines) of radius R are shown. The mem- 
branes are tensionless and the architecture 
/ of the amphiphiles is indicated in the key. 
The instability of the fusion pores at small 
radius is indicated by arrows. For / = 0.3, 
0.35, and 0.4, the stalk-like structure con- 
verts into a pore when it expands to a ra- 
dius R ~ 2ARg at which the free energies 
of stalk-like structure and pore are equal. 
For the system composed of amphiphiles of 
/ = 0.25, however, the stalk-like structure 
converts at i? ~ 2Rg into an inverted micel- 
lar intermediate, IMI, whose free energy is 
shown by the dotted line. The fusion pore 
is unstable to this IMI intermediate when 
its radius decreases to i? ~ 'iARg. Thus the 
IMI is the most stable structure under these 
conditions. From Ref. [ I311j . 

large radii the free energy linearly increases with i?, the slope being proportional to the effective 
line tension of the pore's rim. For very small radii, however, the free energy also increases 
as we decrease the size of the fusion pore because the head groups begin to repel each other 
across the pore. Thus, the SCF calculations predict a barrier for the closing of a fusion pore. 
This prediction qualitatively agrees with the experimental observation of flickering fusion pores 
at very low tension. In this case, fusion pores once formed do not expand due to the lack of 
tension, but they do not close either because of the barrier just mentioned. Thus they remain 
in a metastable state and their radius fluctuates around its preferred value. Experimentally this 
is detected by a flickering of the current between the two fusion vesicles [1323 , 424, 425|. The 
region where flickering pores are expected is also indicated in the fusion diagram (cf. Fig. 121(1 . 

3.4.2. Barriers along the stalk-hole path observed in the simulations 

The fusion path observed in the simulation differs from the classical hypothesis by the oc- 
currence of non-axially symmetric intermediates and the formation of holes in the individual 
bilayers. 

First, we discuss the free energy of an isolated hole in a single membrane. Hole formation 
has been investigated by previous simulations [ 1^ [TT)^ [TT?l IT?^ 1^ II?! 1^ 1^ 1^ H!m 
14321 11751 I433j . These studies have focused on the fluctuations of the pores, their shape and 
detailed structure, e.g., the bulging at the rim (see Fig. [221 middle panel). The results indicate 
that hole formation is well described by classical nucleation theory [ 4511 . Fig. 1231 presents 
the free energy of a hole as a function of its radius R measured in units of the radius of gyration 
Rg of the amphiphilic molecules. While previous SCF calculations [ I307j conceived holes as 
weakly segregated equilibrium structures (the hexagonally perforated phase) of tense bilayers 
the holes that form at large incompatibility must be stabilized by an external constraint (see 
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Figure 23. (a) Free energy of a hole in an isolated bilayer as a function of R/Rg at zero tension 
for various amphiphile architectures, /. From top to bottom the values of / are 0.29, 0.31, 0.33, 
and 0.35. (b) Same as above, but at fixed / = 0.35 and various tensions 7/7int- From top to 
bottom, 7/7mt varies from 0.0 to 0.6 in increments of 0.1. From Ref. [ I312j . 
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Sec. 12. 3() and correspond to unstable structures, in agreement with simulation and experiment. 
For a tensionless membrane, 7 = 0, the free energy linearly increases with its radius, and from 
the asymptotic behavior one can identify the line tension, A of the hole's rim. Upon increasing 
the membrane tension, the free energy curves adopt a parabolic shape expected from classical 
nucleation theory, -Fhoie(-R) = 2ttRX — ttR'^^ [ 434] I435j . Note that the free energy barrier, 
AF* = 7rA^/7, is on the order of a few tens of the thermal energy scale, ksT, for the parameters 
utilized in the MC simulations of the bond fluctuation model. Although tense membranes are 
metastable and eventually will rupture, the high barrier makes the homogeneous nucleation of 
a hole in a bilayer an unlikely event on the time scale on which fusion occurs. Thus, isolated 
bilayer membranes are stable (cf. also Fig. ^ 

In contrast to the above calculation of the free energy of a hole in a bilayer, the direct 
calculation of the free energy of the stalk-hole intermediate is much more difficult. The reason 
is that whereas the former structure is axially symmetric, the latter is not. As noted, the axially 
symmetry is broken explicitly. This means that one has to calculate the free energy of an 
intrinsically three-dimensional structure. This is computationally very demanding. Moreover, 
the choice of a suitable reaction coordinate is less obvious for the complex, non-axially-symmetric 
intermediates of the stalk-hole mechanism. For these practical reasons we estimate the free 
energy of the intermediates observed in the simulations by patching together axially symmetric 
configurations calculated within SCF theory. Two structures, the elongated stalk and the stalk- 
hole complex, are considered. 

The geometry of the elongated stalk, which is observed before a hole forms next to it, is 
characterized by the arc of a circle of radius, R, and length lirRa, with < a < 1. We 
decompose the free energy of this structure into the contribution of its two end caps, which each 
resemble half of a metastable, circular stalk, and a fraction, a, of a ring-shaped stalk, or inverted 
micellar intermediate (IMI) (see Fig]24|). 



By constructing the free energy via the simple addition of the free energy of axially symmetric 
structures and not allowing for additional optimization of the shape, we will overestimate the 
free energy of the intermediate. 

After a hole forms next to the stalk, the geometry of the stalk-hole complex can be character- 
ized by the radius of the hole, R — 5, where 6 is radius of the metastable stalk, and the fraction a 
of the hole's rim that is covered by the stalk. The geometry of the latter configuration resembles 
a hemifusion diaphragm. Thus we can approximate the free energy of the stalk-hole complex by 



where we neglect the free energy costs of the two end points of the stalk. The free energy of 
these two saddle-shaped point defects presumably is small. 

These free energy estimates allow us to explore the free energy landscape, F{R, a) = min [F^s, Kjh] 
of the fusion process as a function of the two reaction coordinates, the radius of the elongated 
stalk, R and the fraction a of the hole's rim covered by the stalk. Here, we additionally have 
assumed that a pore will form with a negligible barrier when the stalk- hole free energy, Fgh, is 
lower than the free energy of an elongated stalk, F^s- The formation of a hole in the vicinity of 
a stalk produces a ridge in the free energy landscape. This is, of course, a simplification, but 
we expect the corresponding free energy barrier to be small. The main effect of this simplifi- 
cation is that the locus of transitions between the extended stalk and the stalk-hole complex 
is a sharp ridge, while in a more complete description it would occur over a range of values 



-^es = -F'stalk + aFiMl{R). 



(70) 



i^sh = aFuM + (1 - a)Fhoie(i? - S) 



(71) 
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Figure 24. (a) Density profile of an IMI - an elongated stalk that closes in itself forming a 
circular structure. The amphiphiles are characterized by / = 0.3. The radius of the IMI, in 
units of the radius of gyration, R/Rg is 3.4. (b) Free energy of an IMI as a function of R/Rg at 
zero tension for various amphiphile architectures, / From Ref. [ I312j . 
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Figure 25. Four free energy landscapes (in units of ksT) of the fusion process, plotted as a 
function of the radius, R (in units of Rg) and circumference fraction a. The architecture of the 
amphiphiles, /, and the value of the tension, j/^int, ai^e given. The dotted line shows a ridge of 
possible transition states, separating two valleys. The region close to the a = line corresponds 
to a barely elongated stalk intermediate (see Eq. ()7U() ). The other valley, close to a = 1 states, 
corresponds to a hole almost completely surrounded by an elongated stalk. The saddle point 
on the ridge, denoted by a white dot, corresponds to the optimal, lowest free energy, transition 
state. From Ref. [ l3T2] . 



Coarse-grained models for biological and synthetic membranes 



63 



over which the free energy of the two competing structures differs by order ksT. Once the hole 
has formed, the free energy decreases as the stalk surrounds the hole (increasing a), and the 
stalk- hole intermediate expands (increasing R). 




1.29 D.31 0.33 0.35 



Figure 26. Plot of a*, which corresponds to the optimal transition state in the stalk- hole 
mechanism, as a function of architecture of the amphiphiles and the tension of the membrane. 
Small values of a* (dark) correspond to leaky fusion events. From Ref. [ I312j . 



The minimum of the free energy along the ridge Fcs{R,ct) = K;h(i?, a) defines the transition 
state characterized by R* and a*. The value, a*, is an important characteristic of the fusion 
intermediate. The larger a* , the larger is the fraction of the hole's rim that is surrounded by the 
stalk when the hole is formed. Thus, larger values of a* correspond to the tighter fusion events 
while small values of a* suggest that fusion is leaky and possibly competes with rupture. The 
value of a* as a function of tension, 7, and amphiphilic architecture, /, is shown in Fig. [201 Small 
values of the tension and small values of /, i.e., large negative curvatures, favor tight fusion. 
Stalks expand more readily because the line tension of the stalk decreases with / (cf. Fig. I23|) : 
the line tension, Xh, of a hole increases for small /; and holes expand less quickly for smaller 
tension. 

In Fig. |27| we compare the free energy of the fusion intermediates along the classical hemi- 
fusion path with that of the stalk-hole mechanism as a function of tension and architecture. 
In agreement with experimental observation, the fusion barrier in both scenarios is the lower 
the more negative is the spontaneous curvature of the monolayers and the larger the membrane 
tension is. For the parameters studied, the stalk-hole mechanism has a slightly smaller free 
energy barrier than that of the classical hemifusion path. The free energy difference, however, 
amounts only to a few ksT indicating that the stalk-hole mechanism is a viable alternative to 
the classical hemifusion path, but that the choice of the pathway will depend on specific details 
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Figure 27. Free energy barriers measured relative to the initial metastable stalk, in units of 
ksT, in (a) the new stalk-hole complex mechanism, and (b) the standard hemifusion mechanism. 
From Ref. [ l3T2] . 
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of the system. This is in agreement with the simulation data. 

4. Conclusion and outlook 

We hope that we have demonstrated that coarse-grained models are a valuable tool for in- 
vestigating universal collective phenomena in bilayer membranes. These models bridge the gap 
between atomistic simulations that are limited to rather small time and length scales, and phe- 
nomenological models that ignore much of the internal structure of the bilayer. They allow for 
direct insights into processes that involve many lipids and take place on time scales of millisec- 
onds and length scales of micrometers. 

With a specific example, the fusion of membranes, we have illustrated what can be learned from 
studying coarse-grained models. The apparent universality of membrane fusion that is observed 
experimentally and the separation of time scales between the local structural relaxation in an 
individual membrane and the fusion process make it a suitable subject of such models. A rather 
simple model already yields direct information about the fusion pathway that cannot easily 
be obtained otherwise, and leads to experimentally verifiable predictions. For high membrane 
tension and asymmetric lipid architecture, the SCF calculations predict transient leakage that 
is correlated in time and space with the fusion event. This is observed in simulation of coarse- 
grained models [ I14()| I16H IM] I414| I415| l4(),Sj that substantially differ in their interactions and 
is corroborated by electrophysiological experiments [ 13881 1590 . Moreover, the coarse-grained 
models resolve the overt contradiction between the stability of isolated membranes or vesicles 
that is necessary for their function and the fusion process that involves hole formation in the 
two apposed membranes. This can be understood as follows. The line tension of holes in the 
membranes is reduced by the presence of stalks. Therefore, the rate of hole formation is dictated 
by heterogeneous nucleation at stalks, which is much greater than the rate of homogeneous 
nucleation in an unperturbed bilayer. The study of fusion also illustrates the fruitful interplay 
between particle-based simulations and field-theoretic techniques. The latter technique provides 
quantitative free energy estimates of the fusion barrier and permits the exploration of a wide 
range of molecular architectures and membrane tensions. It reveals that fusion is strongly 
suppressed outside a narrow range of monolayer curvatures (see "phase diagram" in Fig. I21|) . 
This result suggests ways to control fusion by tuning the stability of stalks and the formation of 
holes. 

The increase of computing resources and algorithmic advances will allow atomistic simula- 
tions to explore phenomena at ever larger time and length scales. The role of coarse-grained 
models, however, will not merely consist in studying problems that are not yet accessible by 
atomistic simulations. Rather the reduction of the complexity of the system will highlight the 
relevant ingredients and allow for a systematic quantitative analysis. Thus, coarse-grained mod- 
els will provide insight into specific mechanisms and principles as well as the degree of universal- 
ity. Coarse-grained models can explore a variety of conformational properties, thermodynamic 
quantities (such as phase behavior, bending elastic constants, tension and areal compressibility) 
and dynamic characteristics (e.g., structural relaxation of the bilayer, dynamics of undulations, 
lateral diffusion of lipids) simultaneously, and the parameters of the model system can be varied 
independently and over a wide range. The wealth of information these models provide concerning 
equilibrium and dynamic properties can be compared both to experiment and to simpler phe- 
nomenological approaches in order to assess the validity of the model description. In this way, 
coarse-grained models identify interesting parameter regimes, test phenomenological concepts, 
and permit systematic exploration of collective phenomena in biological membranes. 

A variety of coarse-grained models has been devised in recent years which differ in the degree 
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of molecular detail, the type of interactions, and the representation of the solvent. In view of 
the richness and complexity of collective phenomena in membrane physics, no single model will 
emerge that captures all aspects. "Systematic" coarse-graining procedures [ IT^ l21Ul [TH [T31 [TBI 
[T7| or tightly coupled multi-scale simulation techniques [ \2W\ I436[ I25H I229[ I437j will provide 
information about the specific interactions that have to be incorporated into coarse-grained 
models to describe the phenomena of interest. 

There are a host of biological problems that can be studied by the kind of coarse-grained 
models we have discussed, and we expect their use in biological physics to increase greatly in the 
near future. The formation of bilayers and vesicles [ [IMl IIIIl ESHl USS IHIl ESI EHU has been 
investigated. The studies of different bilayer phases [ 11871 11801 12U61 11911 1175j . of the kinetics 
of phase transition [ I2()7| I2U8| , of self-assembly [ I17H |438' 4391 1411j , and of phase separation 
in membrane consisting of different components [ IHUl 1114. .4411 IT^ IH^ Bill IH^ BIHl 
I447L I233j have already yielded important information about these fascinating systems. Other 
important questions that can be tackled by coarse-grained models include the interaction between 
membranes, between membranes and solid substrates, and the interplay between membrane 
collective phenomena and peptides, proteins and polymers [ 14481 IUHl I45UI 14521 11831 1453] . 
or active components (e.g., ion pumps [ 1454| Il55l 1156, 457 ^ I183j . or fusion proteins [ 14041 1158] ) . 

We close by considering a simple, but important, application to an area which has received 
much study recently; the possibility of phase separation in the plasma membrane. The basic 
idea is that the plasma membrane enclosing the cell, rather than being a homogeneous mixture 
of lipids and cholesterol, may in fact be rather inhomogeneous, with regions of saturated lipids, 
like sphingomyelin, and cholesterol, aggregated and floating, like "rafts", in a sea of unsaturated 
lipids, like most phosphatidylcholines. The field grew quickly when it was observed that some 
signaling proteins preferred the raft environment. As a consequence, these proteins were not 
distributed randomly throughout the plasma membrane, but were concentrated in the rafts and 
could therefore function more efficiently. 

This seems to be a system which could be studied by coarse-grained models [ 1459j . One object 
would be to capture the relevant properties of cholesterol which appears to order the saturated 
lipids in its vicinity [ 1460j . With a suitable model, one can then calculate phase diagrams and 
compare them to the several experimental ones. 

As interesting as this would be, and it would be very interesting indeed, one could then try 
to understand the observation that rafts act as a nucleation site for the fusion process. One sees 
how "...way leads onto way" and to the study of ever more complex phenomena, those which are 
characteristic of biological systems. 
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